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Abstract 

Motivated  by  the  demands  of  simulating  flapping  wings  of  Micro  Air  Vehicles,  novel 
numerical  methods  were  developed  and  evaluated  for  the  dynamic  simulation  of  mem¬ 
branes. 

For  linear  membranes,  a  mixed-form  time-continuous  Galerkin  method  was  em¬ 
ployed  using  trilinear  space-time  elements.  Rather  than  time-marching,  the  entire 
space-time  domain  was  discretized  and  solved  simultaneously.  Second-order  rates 
of  convergence  in  both  space  and  time  were  observed  in  numerical  studies.  Slight 
high-frequency  noise  was  filtered  during  post-processing. 

For  geometrically  nonlinear  membranes,  the  model  incorporated  two  new  schemes 
that  were  independently  developed  and  evaluated.  Time  marching  was  performed  us¬ 
ing  quintic  Hermite  polynomials  uniquely  determined  by  end-point  jerk  constraints. 
The  single-step,  implicit  scheme  was  significantly  more  accurate  than  the  most  com¬ 
mon  Newmark  schemes.  For  a  simple  harmonic  oscillator,  the  scheme  was  found  to  be 
symplectic,  frequency-preserving,  and  conditionally  stable.  Time  step  size  was  lim¬ 
ited  by  accuracy  requirements  rather  than  stability.  The  spatial  discretization  scheme 
employed  a  staggered  grid,  grouping  of  nonlinear  terms,  and  polygon  shape  functions 
in  a  strong-form  point  collocation  formulation.  The  observed  rate  of  convergence  was 
two  for  both  displacement  and  strain.  Validation  against  existing  experimental  data 
showed  the  method  to  be  accurate  until  hyperelastic  effects  dominate. 
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NOVEL  DISCRETIZATION  SCHEMES  FOR  THE  NUMERICAL  SIMULATION 

OF  MEMBRANE  DYNAMICS 

I.  Introduction 

The  advent  of  the  Micro  Air  Vehicle  (MAV)  has  forced  the  scientific  community 
to  focus  on  a  familiar  but  extremely  complex  mechanism:  the  flapping  wing.  Well 
suited  for  low-speed  or  even  hovering  flight,  the  flapping  wing  provides  thrust,  lift, 
and  control.  It  is  typically  highly  flexible  and  non-uniform  materially  and  geomet¬ 
rically.  Besides  undergoing  large  and  rapid  rigid-body  motion,  it  also  experiences 
large  deformations.  The  true  aeroelastic  system  is  therefore  highly  nonlinear  with 
coupled  structural  and  aerodynamic  phenomena  [27,  45].  The  aerodynamics  alone 
have  proved  more  difficult  to  model  than  those  of  a  fixed  wing  [93] . 

Micro  Air  Vehicles  often  utilize  flapping  wings  with  membranes  as  the  lifting  sur¬ 
faces.  Structurally,  a  common  flapping  wing  design  consists  of  a  main  spar,  stiffeners, 
and  a  thin  flexible  membrane  to  serve  as  the  lifting  surface.  Fabricated  MAVs  tend  to 
employ  much  simpler  variations  on  the  theme.  The  spar  and  stiffeners  (often  called 
battens)  have  bending  and  torsional  stiffness  and  can  be  considered  to  behave  as 
beam  elements,  leading  to  the  current  nomenclature  “beam-membrane  structure.” 
They  preserve  the  form  and  provide  the  forcing  mechanism,  while  the  membrane 
provides  the  lightweight  lifting  surface. 

As  interest  in  flapping  membrane-based  wings  intensifies,  requirements  for  numer¬ 
ical  simulation  will  grow.  Lessons  learned  from  two-dimensional  airfoil  studies  will 
naturally  progress  to  fully  three-dimensional  subjects  for  validation  and  application. 
As  the  computational  expense  grows  accordingly,  demand  will  also  rise  for  novel  and 
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efficient  discretization  schemes  that  retain  or  even  improve  sufficient  solution  accu¬ 
racy.  This  forecasted  need  motivated  the  present  research,  drove  its  objectives,  and 
precisely  defined  its  scope. 

Background  and  Motivation 

Compared  to  a  rigid  airfoil,  a  flexible  membrane’s  passive  reactions  to  aerodynamic 
loads  can  improve  gust  response  [87]  and  stall  characteristics  [119].  Wind  tunnel  ex¬ 
periments  of  a  membrane  stretched  over  a  fixed  frame  have  revealed  the  significance 
of  at  least  five  vibration  modes,  with  strong  interactions  between  the  membrane  os¬ 
cillations  and  vortex  shedding  [113].  Gordnier  [45]  performed  high-fidelity  numerical 
simulations  of  a  similar  configuration  and  observed  both  standing  and  travelling  wave 
responses  in  the  membrane  depending  upon  the  flow  regime,  leading  him  to  conclude 
that  an  advanced  multidisciplinary  approach  was  necessary  to  fully  grasp  the  com¬ 
plicated  system. 

Both  air  loads  and  inertial  forces  contribute  to  structural  deformation.  To  exam¬ 
ine  their  relative  significance,  Yin  and  Luo  [143]  defined  the  mass  ratio  m* ,  which  in 
physical  terms  is  the  ratio  of  inertial  force  of  the  wing  and  the  aerodynamic  pressure. 
They  referred  to  an  earlier  study  of  the  hawkmoth  [27] ,  where  m*  =  5  and  deflections 
were  similar  in  a  vacuum  and  in  air,  indicating  that  inertial  forces  dominate.  In  con¬ 
trast,  dragonfly  wings  are  much  lighter  at  m*  =  1  and  aerodynamic  forces  dominate. 
Another  study  of  the  hawkmoth  [121]  showed  a  23%  difference  in  modal  frequencies 
between  air  and  vacuum  experiments,  leading  them  to  conclude  some  aerodynamic 
coupling  was  present.  Hawkmoth  wings  have  been  found  to  deform  primarily  due  to 
inertial  forces,  and  damping  could  approximate  the  aerodynamic  effects  [27].  Clearly 
the  structural  dynamics  are  a  critical  part  of  this  complex  system,  and  must  be  suit¬ 
ably  modeled  to  achieve  a  useful  outcome. 
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In  nature,  many  insects  [139]  and  bats  [128]  have  beam-membrane  wings.  Because 
they  have  mastered  flight  through  evolutionary  development,  these  living  MAVs  have 
been  the  subject  of  many  studies.  Attempts  to  capitalize  on  their  success  despite 
our  limited  knowledge  have  also  led  to  many  “bio-inspired”  designs.  One  of  the 
most  common  subjects  is  the  hawkmoth  ( Manduca  sexta).  Its  relatively  large  size  is 
helpful  because  it  enables  detailed  experimental  studies  [138],  and  evokes  a  suitable 
size  for  fabricated,  missionized  vehicles  with  payloads.  Its  flapping  mechanism  - 
forcing  only  at  the  shoulder  with  passive  control  through  venation  variations  -  avoids 
the  complex  control  mechanisms  of  actively-deformed  wings  like  those  in  birds  and 
bats.  Conventional  structural  analysis  can  therefore  be  performed  on  the  wings.  The 
prevalence  in  the  literature  of  the  hawkmoth  and  the  models  it  inspired  motivates 
this  study  to  focus  on  membrane  systems  employed  by  this  class  of  MAVs. 

Research  Goals 

As  the  demand  for  higher-fidelity  aeroelastic  simulations  of  membrane-based  wings 
continues  to  grow,  so  will  the  need  for  practical  structural  models  that  are  capable 
of  capturing  the  vitally  important  membrane  dynamics.  In  anticipation  of  this  need, 
the  purpose  of  this  study  was  to  devise  and  evaluate  numerical  schemes  to  accurately 
simulate  the  dynamics  of  a  membrane  like  those  employed  in  typical  Micro  Air  Ve¬ 
hicle  wings.  Only  the  membrane  component  was  addressed  in  this  study  because,  as 
mentioned  previously,  the  complexity  of  an  entire  flapping  wing  system  requires  an  in¬ 
terdisciplinary  approach.  The  numerical  schemes  developed  here  will  be  demonstrated 
using  fixed  or  prescribed  frames  only.  Having  thoroughly  verified  and  validated  the 
results,  future  incorporation  of  the  schemes  into  flapping  beam-membrane  structures 
is  a  natural  progression. 

The  scope  and  requirements  for  this  study  will  now  be  delineated.  The  hyperbolic 
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governing  partial  differential  equations  (PDEs)  for  a  membrane  are  second  order  in 
space  and  time  wave  equations).  All  proposed  schemes  were  to  be  capable  of 
simulating  the  transient  response  of  a  membrane,  including  standing  and  travelling 
waves.  Modal  analysis  was  therefore  expressly  excluded  -  if  needed,  modal  frequencies 
could  be  extracted  from  the  time  series  information  of  the  transient  response.  Accu¬ 
racy  was  to  meet  or  exceed  that  of  commonly-used  numerical  computation  methods. 

Capability  requirements  and  scope  of  suitability  for  the  new  numerical  schemes 
were  derived  from  the  characteristics  of  the  physical  system  under  consideration:  a 
fabricated  (as  opposed  to  biological)  structure  consisting  of  a  thin,  isotropic  membrane 
stretched  across  a  rigid  frame.  Selecting  this  system  enabled  the  design  to  focus  on 
a  narrow  but  useful  design  space.  The  model  was  therefore  customized  for  systems 
abiding  by  the  following  assumptions: 

•  Geometry 

—  Large  deformation/small  strain  (at  least  geometrically  nonlinear) 

—  Variable  geometry  (not  restricted  to  circle,  square,  etc.) 

—  Membranes  must  be  planar  when  at  rest 

•  Membrane  properties 

—  Isotropic 

—  Buckling  and  wrinkling  were  neglected 

—  Linear  clastic  material  model 

—  Poisson  ratio  u:  0  to  0.5  (hawkmoth  approximation  v  =  0.49  [27]) 

—  Membrane  at  rest  may  be  slack  or  prestressed 

—  Clamped  at  all  boundaries  (Dirichlet  boundary  conditions  only) 
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Forcing  functions 


—  All  external  and  body  forces  are  distributed  and  smooth  (discontinuous 
point  forces  not  expected) 

—  Directional  forces  (tractions  or  body  forces  like  gravity) 

—  Follower  forces 

—  Tightly-coupled  aerodynamic  forces  that  respond  to  the  membrane  shape 

—  Time-dependent  forces 

•  Frame 

—  Prescribed  motion  only.  Incorporation  of  a  coupled  beam  model  for  the 
rigid  frames  was  beyond  the  scope  of  this  study. 

The  performance  of  the  proposed  techniques  was  assessed  through  a  rigorous  pro¬ 
gression  of  verification  and  validation.  As  precisely  defined  in  Ref.  [112],  the  veri¬ 
fication  process  ensures  that  a  model  is  solving  the  chosen  equations  correctly  and 
consistently.  After  code  has  been  verified,  validation  addresses  whether  the  proper 
equations  were  selected  for  approximating  the  physical  system  under  consideration. 
Experimental  data  therefore  plays  a  key  role  in  validation. 

A  wide  variety  of  methods  for  building  a  simulation  that  meets  the  listed  require¬ 
ments  are  available,  and  many  of  them  will  be  discussed  in  the  literature  review  in 
Chapter  Two.  For  two-dimensional  membranes  of  arbitrary  geometry,  the  Galerkin 
finite  element  method  (FEM)  is  the  predominant  approach.  The  most  commonly  as¬ 
sociated  time  integration  scheme  is  the  Newmark  family  of  methods,  favored  mostly 
for  the  unconditional  stability  of  its  implicit  average  acceleration  scheme.  Because 
of  the  Galerkin  formulation’s  numerical  integration  procedures,  the  computational 
expense  of  this  approach  for  nonlinear  dynamic  problems  is  significant.  In  addition, 
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for  highly  nonlinear  problems,  the  stability  benefits  of  the  implicit  Newmark  methods 
are  degraded. 

The  narrowed  design  space  for  this  study,  however,  provided  the  opportunity 
to  develop  new  methods  tailored  to  the  requirements.  Hence,  novel  schemes  for 
both  time  and  space  discretization  were  created,  evaluated,  and  incorporated  into 
a  final  membrane  model.  In  particular,  a  time-marching  scheme  based  on  Hcrmite 
polynomial  interpolation  was  developed  that  is  significantly  more  accurate  than  the 
second-order  Newmark  methods,  yet  it  is  stable  at  time  step  sizes  appropriate  for 
capturing  the  relevant  dynamics.  Likewise,  a  spatial  discretization  scheme  based 
on  the  point  collocation  method  was  devised  that  handles  the  geometric  nonlinear 
membrane  behavior  while  avoiding  the  requirement  for  numerical  integration.  When 
put  together,  these  two  schemes  provided  a  practical,  robust  model  for  the  dynamic 
simulation  of  geometric  membranes. 

Organization 

Thus  far  the  motivating  systems  and  the  scope  of  the  study  have  been  outlined. 
The  remainder  of  the  dissertation  will  detail  the  individual  components  of  the  model 
and  finally  evaluate  the  complete  membrane  model.  The  literature  review  will  be 
conducted  in  Chapter  Two  to  frame  the  new  methods  in  the  context  of  conventional 
approaches.  Chapter  Three  details  the  development  and  evaluation  of  the  Simulta¬ 
neous  Time-Continuous  Galerkin  (STCG)  method  for  a  linear  membrane.  Chapter 
Four  progresses  to  fully  nonlinear  membranes.  The  jerk-based  constraint  formulation 
for  Hermite  time  interpolation  is  proposed  and  evaluated,  followed  by  a  similar  devel¬ 
opment  of  the  staggered-grid  point  collocation  model.  Having  rigorously  examined 
these  two  schemes  individually,  they  are  finally  incorporated  into  the  final  membrane 
model  for  analysis.  Lastly,  in  Chapter  Five,  results  will  be  summarized  and  discussed 
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from  the  perspective  of  the  study  objectives.  Recommendations  for  further  study  will 
also  be  offered. 
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II.  Literature  Review 


What  exactly  is  a  membrane?  Jenkins  provides  excellent  definitions  depending 
on  one’s  perspective:  a  membrane  is  a  structure  that  does  not  resist  bending,  while 
a  membrane  model  is  “an  idealized  model  of  a  plate  or  a  shell  structure,  wherein  the 
in-plane  response  dominates  away  from  domain  and  load  boundaries.”  [63].  With  its 
total  lack  of  bending  stiffness,  a  membrane  may  be  thought  of  as  a  two-dimensional 
version  of  a  string  or  cable  [95].  When  subjected  to  an  external  load,  the  only  restoring 
force  is  from  the  tension  in  the  plane  of  the  membrane.  Membranes  are  inherently 
nonlinear  structures  that  may  undergo  large  rigid  body  displacements  and  rotations, 
and  often  large  strains  as  well.  Thorough  discussions  of  membrane  mechanics  and 
challenges  may  be  found  in  [64],  [63],  [96],  and  [95]. 

The  prototypical  problem  for  static  membrane  study  is  the  circular  membrane 
with  fixed  circumference  that  is  inflated  under  constant  pressure.  Commonly  called 
the  Hencky  problem  [51],  this  configuration  facilitates  experimentation  and  analytical 
solutions.  Finite  element  (FE)  models  have  been  applied  to  the  Hencky  problem  that 
use  detailed  geometric  surface  descriptions  with  Jaumann  strains  and  stresses  [96]. 
Another  approach  is  to  use  simple  linear-elastic  elements  applicable  to  the  large  dis¬ 
placement/small  strain  regime  [97]. 

Despite  the  intrinsic  nonlinearity  of  true  membranes,  simplifying  assumptions  may 
be  appropriate  for  some  applications.  Based  on  these  possible  assumptions,  a  mem¬ 
brane  model  can  be  classified  as  “linear,”  “small  strain/finite  rotation,”  or  “fully 
nonlinear”  based  on  how  (or  whether)  it  handles  the  geometric  and  material  non- 
linearities  [64],  Stanford  et  al.  [125]  referred  to  these  three  categories  as  “low-”, 
“medium-”,  and  “high-fidelity.” 


•  Linear  (Low  fidelity) .  This  model  is  the  simplest,  but  carries  the  most  stringent 
assumptions.  The  material  is  linear  clastic.  Initial  membrane  internal  forces 
(pre-tension)  are  significantly  greater  than  any  change  in  internal  forces  caused 
by  deformations.  External  loads  are  not  affected  by  the  membrane’s  shape. 
Only  transverse  displacements,  w,  are  permitted,  and  the  displacements  are 
small.  The  linear  model  may  be  appropriate  only  in  cases  with  small  pressures, 
small  displacements,  and  large  initial  tensions  [126,  125].  The  linear  model  is 
invalid  (unbounded)  for  a  slack  membrane  [125].  Further,  it  has  been  found 
that  a  circular  membrane  that  is  planar  when  undeformed  cannot  be  linearized, 
no  matter  how  small  the  displacements  [88]. 

•  Small  strain/finite  rotation  (Medium  fidelity).  Also  called  a  “geometrically 
nonlinear”  model.  In  this  case,  the  membrane’s  resistance  to  an  increase  in 
transverse  load  (a  force  normal  to  the  plane  of  the  membrane)  depends  on  the 
current  shape  of  the  membrane.  For  example,  if  the  membrane  is  flat,  the 
membrane  offers  no  resistance  when  a  load  is  applied.  On  the  other  hand,  if 
the  membrane  is  deformed  out-of-plane,  the  surface  is  sloped  such  that  tension 
forces  are  able  to  resist  further  deformation  when  the  load  is  incremented.  This 
phenomenon  occurs  even  though  the  material  is  considered  to  be  linear  elastic. 

•  Fully  nonlinear  (High  fidelity).  In  a  fully  nonlinear  model,  both  material  and 
geometric  nonlinearities  are  accommodated.  Often,  the  material  is  considered  to 
be  hyperelastic  -  the  strain  energy  depends  on  both  the  initial  and  the  deformed 
state,  but  is  not  dependent  on  the  path  between  those  states  [15]. 

Geometrically  nonlinear  models  have  been  shown  to  be  suitable  until  hyperelastic 
effects  begin  to  dominate.  For  an  initially  flat,  circular,  rubber  membrane  with  a 
clamped  circumference,  this  transition  occurs  at  a  center  deflection  of  approximately 
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25%  of  the  membrane’s  radius  [106].  Since  this  degree  of  deformation  was  consid¬ 
ered  unlikely  for  typical  Micro  Air  Vehicle  wings,  the  fully  nonlinear  model  was  not 
considered  in  this  study. 

With  the  fundamental  characteristics  of  membranes  in  mind,  the  rest  of  this  chap¬ 
ter  reviews  numerous  approaches  for  modeling  a  membrane.  Broadly  speaking,  they 
can  be  categorized  as  analytical  or  numerical.  Methods  for  handling  the  temporal 
dimension  in  a  numerical  scheme  are  also  discussed. 


Analytical  Models 


Analytical  solutions  for  membranes  are  restricted  to  special  subcases  where  forcing 
functions,  initial  conditions,  and  boundary  conditions  are  tractable.  As  a  result, 
configurations  are  typically  rectangular  for  Cartesian  coordinates,  or  circular  for  polar 
coordinates  or  to  capitalize  on  radial  symmetry. 

The  linear  membrane  model  is  most  amenable  to  analytical  treatment,  and  is 
therefore  most  prevalent  in  the  historical  and  canonical  literature  [91,  136,  46,  108, 
133].  The  governing  equations  may  be  obtained  by  either  a  variational  approach 
using  energy  formulations  and  Hamilton’s  Principle,  or  a  Newtonian  force  equilibrium 
approach.  The  membrane’s  surface  is  defined  by  the  displacement  w(x,y,t)  of  a 
particle  ( x ,  y)  in  the  z  direction.  Letting  p  be  the  membrane’s  constant  density  (mass 
per  unit  area),  P  be  the  constant  tension  (force  per  unit  length),  and  f(x,y,t )  be 
the  external  pressure  in  the  same  direction  as  w  (force  per  unit  area),  the  well-known 
second-order  PDE  is 

Cp'll) 

P-QP  =  PV2™  +  /  (1) 

with  the  Laplacian  operator  defined  as 


82  d2 
dx2  d y2 


(2) 
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The  wave  speed  c,  where  c  =  \JP/p,  is  variously  called  the  the  characteristic  velocity, 
the  speed  of  sound  for  the  material  [28],  or  the  solution  propagation  speed  [58].  For 
the  static  case,  the  Poisson  equation  is  recovered.  If  shear  forces  are  included,  the 
internal  forces  are  represented  by  Nx,  Ny,  and  Nxy,  and  the  static  version  of  Eq.  1 
may  be  presented  as  [125] 


AT  d2w 

NX^T  +  Nr 


dx 2 


xy 


d2w  AT  d2w 

+  Ny^rr^  +  /  =  0 


dxdy 


dy 2 


(3) 


Exploiting  the  radial  symmetry  of  circular  membranes  permits  a  two-dimensional  de¬ 
velopment  to  solve  a  three-dimensional  problem.  The  complementary  energy  principle 
has  been  used  to  derive  a  fully  nonlinear  model  for  predicting  large  deformations  [88] . 
Likewise,  equilibrium  principles  have  been  employed  with  assumed  functions  for  the 
dependent  variables;  the  solution  is  obtained  after  solving  for  the  constant  coefficients 
multiplying  the  assumed  functions  [35]. 

The  perturbation  method  has  also  been  effective  for  membrane  analysis.  A  per¬ 
turbation  method  of  the  pressure  term  has  been  used  to  analyze  a  small  strain/finite 
rotation  circular  membrane.  The  numerical  solution  was  obtained  through  a  Taylor 
expansion  of  the  analytical  solutions  to  differential  equations  [127].  A  perturbation 
of  the  external  force  to  the  same  problem  has  also  yielded  accurate  results  [32],  For 
the  dynamic  case,  the  vibrations  of  a  large-displacement  nonlinear  square  membrane 
were  accurately  modeled  by  using  the  perturbation  method  to  derive  simplified  ap¬ 
proximate  governing  equations  [25]. 


Numerical  Models 

The  simplifying  assumptions  required  by  analytical  approaches  are  quickly  vio¬ 
lated  when  more  complex  systems  are  considered.  The  inclusion  of  arbitrary  domain 


11 


geometry  or  interactive  external  forces,  for  example,  require  a  numerical  approach 
to  approximate  the  solution.  The  three  primary  approaches  are  finite  differences 
(FD),  finite  volume  (FV),  and  finite  element  analysis  (FEA)  [74],  Although  all  three 
approaches  can  be  universally  applied  to  approximate  the  solutions  to  differential 
equations,  there  are  strengths  and  weaknesses  associated  with  each.  As  a  result,  the 
first  two  approaches  are  generally  employed  by  the  fluid  dynamics  and  heat  transfer 
communities,  and  FEA  is  heavily  favored  by  structural  analysts.  For  this  reason, 
this  dissertation  will  focus  primarily  on  FEA.  However,  certain  applications  of  one 
method  can  lead  to  similarities  (if  not  overlap)  with  another,  so  we  should  resist  the 
tendency  to  pigeonhole  methods  to  one  application  or  another. 

The  computational  expense  of  a  particular  model  is  affected  largely  by  the  phe¬ 
nomena  the  analyst  wishes  to  include.  Thus,  one  should  strive  to  use  the  most  efficient 
path  that  produces  answers  of  sufficient  accuracy.  In  one  extreme,  rigid  plates  have 
been  used  for  kinematic  optimization,  where  the  minimizing  computational  effort  was 
paramount  [13].  This  is  the  exceptional  case,  however  -  estimates  of  aerodynamic 
power,  work,  thrust  are  commonly  sought  and  require  a  non-rigid  body  to  capture 
the  physics.  The  numerical  simulation  of  a  flexible,  dynamic,  nonlinear  structure  is 
a  challenging  endeavour,  and  following  sections  will  review  available  methods. 

Classical  Mechanics  and  Hamiltonian  Systems. 

Classical  mechanics  form  the  mathematical  basis  for  the  variational  development 
of  finite  element  models  and  spring-mass  models.  A  brief  review  will  be  presented  here 
as  a  foundation  for  later  developments.  As  opposed  to  the  force-oriented  Newtonian 
mechanics,  classical  mechanics  apply  variational  calculus  to  scalar  energy  functionals 
to  determine  the  optimal  solution  of  a  system.  The  classical  Lagrangian  formulation 
is  well  known  [42,  78,  44],  and  leads  to  Hamilton’s  canonical  equations.  More  recently, 


12 


direct  solutions  have  been  obtained  using  Hamilton’s  Weak  Principle,  which  is  derived 
from  Hamilton’s  Law  of  Varying  Action  (HLVA)  [80].  Both  of  these  approaches 
have  been  utilized  to  derive  numerical  methods  for  approximating  the  behavior  of 
Hamiltonian  systems. 

The  Lagrangian  L  of  a  mechanical  system  is  defined  as  the  difference  between  the 
system’s  kinetic  energy  T  and  potential  energy  U.  Let  q  be  the  vector  of  generalized 
coordinates  according  to  the  Cartesian  position  X  and  time  t.  Then 

L(t,q,q)  =  T(q)  —  U(q)  (4) 


The  physical  path  of  the  system  from  time  f0  to  time  t\  satisfies  the  Extended  Hamil¬ 
ton’s  Principle  [118,  108], 


5  f  Ldt  + 
Jt0 


SWdt  =  0 


't0 


(5) 


To  transform  the  formulation  from  Lagrangian  to  Hamiltonian,  we  next  define  the 
generalized  momenta  p  of  the  generalized  coordinates: 


Pi  = 


d  L 
dq. 


(6) 


The  Hamiltonian  H  is  then  defined  as 


H  =  Y1  ~  L  (7) 

3= 1 

and  with  our  formulation  of  T  it  can  be  proved  that  H(q,p)  =  T(p)  +  U(q).  Substi- 
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tution  into  the  functional  /  gives  the  equation  in  terms  of  the  Hamiltonian: 


Finding  the  extremum  of  this  functional,  one  obtains  the  Hamilton  equations  of  mo¬ 
tion  (also  called  the  “canonical  equations  of  Hamilton”  [78]). 


H{q,p)  __  H(q,p) 
)  Pi 


(9) 


For  direct  numerical  application,  the  configuration  state  vector  r)  is  defined  to  contain 
the  coordinates  q  and  momenta  p, 


V  = 


M 

M 


(10) 


leaving  H  as  a  function  only  of  rj.  The  partial  derivatives  of  H  with  respect  to  rj  are 


dH/dqi 


dH(rj) 

dq 


I  dHj  dqN  I 
I  dH/dpi  f 


(11) 


dH/dpN 


The  canonical  equations  can  then  be  put  into  matrix  form,  which  is  more  amenable 
to  numerical  implementation. 


0  I  3H(p) 

-I  0  dr> 


(12) 
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The  block  matrix  is  often  labeled  J  (not  to  be  confused  with  a  Jacobian  matrix), 
resulting  in  the  form  [44] 


f]  =  J 


dH(rj) 
dr i 


(13) 


With  simple  difference  equations,  this  form  can  easily  be  turned  into  a  midpoint  value 
integration  scheme  [22,  80], 


m  -  >)o  _  3dHQ^) 
At  dr] 


(14) 


The  second  approach,  direct  solution  based  on  Hamilton’s  Weak  Principle  (HWP), 
retains  boundary  terms  in  Hamilton’s  Law  of  Varying  Action  that  are  neglected  in 
Hamilton’s  Principle  [7].  By  honoring  these  terms,  the  energy  balance  through  a 
time  step  is  preserved.  The  weak  form  captures  all  of  the  boundary  conditions  in  a 
single  functional  and  reduces  the  interpolation  order  requirements.  Constant  shape 
functions  may  be  permissible,  leading  naturally  to  discontinuous  Galerkin  formula¬ 
tions  where  the  integration  can  be  performed  by  inspection  [4,  80].  In  fact,  weak 
Galerkin  formulations  and  HWP  result  in  equivalent  schemes  [18].  Representing  non¬ 
conservative  generalized  forces  as  Q ,  Hamilton’s  Weak  Principle  as  developed  in  [80] 
is 

/  n  n  n  \  n  n 

Y  pMj  -  Y  <iisPj  ~6H+Y  QMj  ) df  ~  51  pMi  K + Y  s A  IS  = 0 

j  \j= i  j= i  j= i  y  j= i  j= i 

(15) 

Statements  of  HWP  may  also  be  found  in  [18,  16,  56,  57]. 


Spring-Mass  Models. 

Spring-mass  models  of  a  membrane  discretize  a  continuum  into  a  system  of  in¬ 
terconnected  particles.  Movement  of  the  particles  determines  the  kinetic  energy  of 
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the  system.  Mass  lumping  at  the  particles  is  intrinsic  (the  consistent  mass  treat¬ 
ment  of  FEM  is  not  an  option)  and  the  rotational  momentum  about  points  between 
the  particles  is  ignored.  However,  for  most  problems  and  with  small  elements,  the 
rotational  inertia  is  insignificant  relative  to  the  translational  inertia  [134],  Massless 
springs  connect  particles  and  house  the  potential  energy.  The  model  intrinsically  han¬ 
dles  geometric  nonlinearities  and  (through  the  spring  model)  material  nonlinearity. 
Systems  can  range  from  a  simple  harmonic  oscillator  to  a  cloth  blanket  draped  over 
a  table.  Spring  mesh  models  offer  the  advantages  of  simplicity  and  computational 
speed  [132],  These  factors  make  them  popular  in  the  held  of  animation,  where  they 
have  found  the  most  use  and  development  [73].  One  must  be  aware,  though,  that 
in  this  milieu  realistic-looking  results  may  be  the  objective  rather  than  physically 
accurate  ones  [72],  Despite  this  caution,  the  held  offers  interesting  alternatives  to 
traditional  techniques  in  solid  mechanics. 

One  study  [73]  deemed  spring-mass  models  to  be  superior  to  traditional  FE  for¬ 
mulations  for  purely  axial  structures  undergoing  large  displacements,  in  particular 
the  hanging  chain  or  net.  Algorithm  details  were  not  presented  other  than  stating 
that  Runge-Kutta  time  integration  was  used.  Demonstrations  of  three-dimensional 
static  problems  showed  the  technique  to  be  promising  but  with  some  limitations  (in 
particular,  computational  expense  and  potential  instability  in  complex  problems). 

According  to  Gelder,  as  of  1998  little  work  had  been  done  to  compare  results 
of  spring  mesh  models  with  those  of  traditional  finite  element  models  [132],  He 
evaluated  the  capability  of  a  spring  mesh  model  to  model  an  isotropic  linearly  elastic 
membrane  undergoing  in-plane  deformations.  All  simulations  were  static  analyses. 
He  compared  the  stiffness  matrices  of  the  spring  mesh  model  and  the  constant  strain 
model  to  prove  the  spring  mesh  cannot  exactly  match  the  constant  strain  model 
when  the  springs  have  the  same  stiffness.  In  particular,  the  spring  model  does  not 
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include  the  Poisson  effects.  Based  on  the  geometry  of  an  element  and  the  desired 
edge  stresses,  he  modified  the  spring  stiffness  coefficients  to  more  closely  match  the 
constant  strain  triangle  model.  Exact  solution  was  still  not  possible. 

Delingette  [31]  rigorously  developed  objects  called  “triangular  biquadratic  springs” 
in  such  a  way  that  the  membrane  strain  energy  was  exactly  that  of  an  equivalent  con¬ 
tinuum  model.  On  an  unstructured  triangular  mesh,  these  springs  effectively  modeled 
static  non-linear  membrane  deformations.  As  contrasted  with  previous  efforts,  this 
study  included  an  angular  stiffness  to  the  springs.  This  mechanism  enabled  inclusion 
of  the  Poisson  effect. 

Finite  Element  Models. 

The  predominant  numerical  method  for  membrane  simulation  is  the  finite  element 
method  (FEM).  For  MAV  applications,  different  approaches  have  included  distinct 
beam  and  shell  elements  [90,  34],  membrane  elements  with  varying  parameters  [86], 
or  plate  elements  with  bending  stiffness  set  to  near  zero  for  membrane  portions  of  the 
structure  [124], 

In  FEA,  membranes  are  commonly  considered  subcases  of  shells  [95]  and  element 
development  follows  that  train  of  thought.  However,  the  nature  of  a  membrane  struc¬ 
ture  (light  weight,  minimal  thickness,  large  displacements  and  rotations,  etc.)  causes 
unique  problems  for  conventional  FEA.  While  most  production  codes  have  membrane 
elements  incorporated,  their  accuracy  may  be  questionable  for  some  cases;  in  partic¬ 
ular,  complex  phenomena  such  as  wrinkling  may  not  be  accurately  predicted  [95].  As 
a  result,  attempts  have  been  made  to  custom  tailor  elements  to  this  problem. 

Pauletti  [97]  refined  a  membrane  element  designed  specifically  for  geometrically 
nonlinear  membranes.  The  original  element  is  called  “Argyris’s  Natural  Membrane 
Finite  Element”  and  is  based  on  the  constant  strain  triangle.  The  model  splits  the 
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element  stiffness  matrix  into  three  components  that  capture  geometric  nonlinearities, 
constitutive  relationships,  and  external  load  effects  separately.  The  membrane  is  as¬ 
sumed  to  be  linearly  elastic.  Comparisons  to  benchmark  cases  demonstrated  good 
performance.  The  linear  membrane  displayed  limited  accuracy  under  large  deforma¬ 
tions.  Leung  [83]  developed  trapezoidal  finite  elements  using  trigonometric  shape 
functions  instead  of  polynomials.  Modal  frequencies  were  obtained  more  accurately 
than  with  the  conventional  polynomial  approach,  and  with  greater  conditioning  and 
stability. 

Clearly  there  is  no  consensus  about  how  to  tackle  membranes  in  FEA,  and  special 
care  must  be  taken  to  capture  the  unique  behavior  of  a  particular  membrane  problem. 

Point  Collocation  and  Group  Finite  Element  Methods. 

The  model  developed  here  will  utilize  two  approaches  which  are  relatively  rare 
in  structural  FEM  applications:  the  point  collocation  method  and  the  group  finite 
element  (FE)  formulation.  Both  approaches  were  selected  for  this  effort  because  they 
tend  to  result  in  simpler  formulations  [26,  146],  potentially  offering  greater  code  flexi¬ 
bility  without  sacrificing  accuracy.  Coincidentally,  both  the  point  collocation  method 
and  group  FE  formulation  have  close  ties  with  the  development  of  computational  fluid 
dynamics  (CFD)  schemes.  Development  of  the  group  FE  formulation  in  the  litera¬ 
ture  centers  around  fluid  mechanics  applications  and  examples  [26,  37,  38],  and  the 
similar  lumping  of  nonlinear  terms  in  the  flux  vector  is  standard  practice  for  finite 
volume  formulations  [55].  Similarly,  the  point  collocation  method  shares  a  history 
with  and  bears  a  resemblance  to  the  finite  difference  method  [146].  The  nonlinear 
wave  equations  that  describe  membrane  dynamics  offer  an  interesting  test  case  for  fur¬ 
ther  examining  how  well  these  CFD-associated  techniques  transfer  into  the  structural 
dynamics  milieu. 
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The  collocation  method  is  commonly  mentioned  in  the  finite  element  literature 
when  listing  the  members  of  the  Method  of  Weighted  Residuals  family.  However,  it 
is  rarely  seen  in  practice.  The  weighting  function  for  each  designated  point  in  the 
domain  is  the  Dirac  delta  function.  By  definition,  the  Dirac  delta  function  equals  zero 
everywhere  except  at  its  associated  point,  and  its  integral  over  the  domain  equals 
one  [146].  The  resulting  system  of  equations  solve  the  PDE  (Partial  Differential 
Equation)  point-wise  rather  than  in  an  integral  sense.  Thus,  there  is  no  need  for 
Gaussian  integration  over  an  element.  Posed  in  the  strong  form,  the  method  requires 
the  interpolation  scheme  to  be  differentiable  to  the  same  order  as  the  PDE  [146,  1], 
but  the  formulation  of  the  system  of  equations  has  been  found  to  be  less  complicated 
than  when  employing  the  Galerkin  method  [1]. 

Also,  while  certainly  applicable  to  a  conventional  FEM  mesh,  the  point  collo¬ 
cation  method  lends  itself  to  a  variety  of  unconventional  discretization  and  inter¬ 
polation  schemes.  Element  shapes  can  expand  beyond  triangles  and  quadrilater¬ 
als  to  n-sided  polygons.  Taking  the  concept  even  further,  “meshless”  [10,  61]  and 
“element-free”  [11,  76]  methods  use  least-squares  fits,  radial  basis  functions,  or  other 
neighboring- node-based  techniques  for  forming  the  system  of  equations.  Using  a 
meshless  point  collocation  method,  [1]  solves  a  wide  variety  of  problems  including 
heat  conduction,  Couette  flow,  and  a  cantilever  beam.  In  the  interest  of  computa¬ 
tional  efficiency  for  the  dynamic  simulation,  the  present  model  utilizes  a  staggered 
background  mesh  so  it  does  not  fall  into  the  meshless  category. 

In  a  group  finite  element  formulation  (also  called  “product  approximation”  [26]), 
aggregated  nonlinear  terms  are  first  computed,  then  interpolated  as  a  single  degree  of 
freedom.  Consider  the  term  puv  where  p,  u,  and  v  are  each  dependent  variables  [38]. 
Rather  than  applying  trial  solutions  <pj  to  each  of  the  variables,  the  group  formulation 
interpolates  their  precomputed  products  as  <t>j  (puv)j-  Significant  computational 
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savings  have  been  observed,  with  the  benefits  increasing  from  higher  dimensionality  or 
order  of  nonlinearity  [37].  For  some  cases,  point-wise  accuracy  may  actually  be  higher 
than  that  of  the  Bubnov-Galerkin  method  [26].  This  observation  hints  at  a  potential 
synergy  in  the  pairing  of  the  group  FE  formulation  with  the  point  collocation  method. 

Time  Integration  Methods 

Numerical  approximation  of  initial  value  problems  (IVP),  also  aptly  called  evolu¬ 
tionary  differential  equations  [3],  has  been  studied  for  centuries,  with  Euler’s  foun¬ 
dational  work  coming  in  the  18th  Century.  Despite  its  rich  history,  the  challenge  of 
balancing  stability  and  accuracy  continues  to  invite  further  consideration,  especially 
as  computational  models  are  applied  to  progressively  more  complex  systems. 

A  primary  source  of  problems  for  evolutionary  equations  is  the  propagation  and 
build-up  of  errors.  Hence,  dissipation  of  unwanted  high-frequency  content  is  desired. 
Ideally  it  is  achieved  through  the  formulation  of  the  predictive  algorithm,  often  em¬ 
ploying  parameters  to  tailor  the  behavior  [54] .  In  some  circumstances,  this  dissipative 
behavior  may  be  necessary  because  noise  can  significantly  affect  extrapolated  predic¬ 
tions  and  result  in  instability  [58].  Spurious  high-frequency  oscillations  have  long 
been  documented  for  Galerkin  solutions  of  the  wave  equation  [92]  and  the  shallow 
water  wave  equations  [5,  39].  The  time-discontinuous  Galerkin  (TDG)  method  offers 
additional  flexibility  through  the  flux  resolution  to  control  dissipation  [53];  upwind- 
ing  techniques  for  the  jump  operators  in  particular  have  proven  effective  in  improving 
stability  [59,  16]. 

Structural  dynamics  formulations  spring  from  the  equilibrium  equation 

Mu  +  Cu  +  Ku  +  F  =  0  (16) 
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where  u  is  the  solution  vector,  M,  C,  and  K  are  the  mass,  damping,  and  stiffness 
matrices  respectively,  and  F  is  the  external  force  vector. 

The  most  common  time  integration  technique  in  FEM  is  direct  integration.  The 
direct  integration  approach  splits  the  domain  into  two  parts:  a  spatial  component 
to  be  addressed  with  a  finite  element  approximation  and  a  time-marching  predictor 
algorithm.  The  FEA  is  a  static  solution  where  for  linear  systems  the  forces  exerted 
by  the  structure  are  a  function  of  the  structure’s  current  configuration.  Given  appro¬ 
priate  initial  conditions,  external  forces,  and  the  internal  forces  from  the  FEA,  the 
time- marching  algorithm  predicts  the  configuration  at  the  end  of  a  time  increment. 
That  solution  then  provides  the  initial  conditions  for  the  next  time  step,  and  so  on. 
Much  effort  has  gone  into  developing  algorithms  that  balance  the  often  contradictory 
demands  placed  upon  them.  Numerous  families  of  algorithms  have  been  developed, 
including  the  single-step  direct  integration  methods  that  are  either  explicit  (such  as 
the  Euler  and  central  difference  methods)  or  implicit  (Newmark,  Wilson,  etc.).  These 
procedures  produce  at  best  second-order  accuracy  [58].  Explicit  methods  are  more 
economical  at  each  step,  but  they  are  conditionally  stable  and  typically  require  more 
steps.  Implicit  methods  require  more  expensive  calculations  but  can  be  formulated 
to  be  unconditionally  stable,  enabling  larger  time  steps. 

Besides  stability  limits,  the  time  step  size  for  nonlinear  systems  is  limited  by 
other  factors.  A  chaotic  system  may  force  strict  accuracy  requirements  to  minimize 
the  propagation  of  numerical  perturbations.  For  implicit  schemes,  iterations  may  not 
converge  if  the  step  size  is  so  large  that  a  sufficiently  good  initial  guess  cannot  be  made. 
These  requirements  may  limit  step  size  to  the  point  that  conditionally  stable  schemes 
become  feasible.  From  these  considerations  arose  the  goal  of  implementing  a  method 
that  offers  higher  accuracy  than  the  predominant  Newmark  methods.  Effective  han¬ 
dling  of  nonlinearities,  conservation  of  energy,  improved  frequency  preservation,  and 
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useful  error  estimation  for  step  size  or  damping  control  were  additional  objectives. 

In  the  context  of  structural  finite  element  analysis,  the  conventional  approach  for 
marching  this  type  of  model  through  time  is  to  use  one  of  the  members  of  the  well- 
documented  Newmark  family  [58,  68,  28].  All  of  the  Newmark  methods  are  single-step 
methods  that  are  distinguished  by  the  values  chosen  for  two  parameters,  f3  and  7. 
The  two  most  commonly  used  methods  are  the  explicit,  conditionally  stable  central 
difference  method  and  the  implicit,  unconditionally  stable  trapezoidal  rule  (also  called 
the  average  acceleration  method).  Both  methods  are  second-order  accurate  in  time. 
Unconditional  stability  is  generally  attractive  for  structural  applications — despite  the 
presence  of  high-frequency  content  or  numerical  noise,  large  time  steps  can  be  taken 
while  preserving  the  accuracy  of  the  relevant  low-frequency  modes.  It  should  be 
noted  that  unconditional  stability  is  no  longer  guaranteed  for  non-linear  problems 
and  the  order  of  the  methods  may  drop  to  one  in  the  presence  of  damping  [28]. 
Despite  the  popularity  of  the  Newmark  methods  in  structural  analysis,  other  time 
integration  techniques  such  as  the  Runge-Kutta  method  and  space-time  finite  element 
formulations  have  been  successfully  applied,  showing  there  is  room  for  alternatives. 

The  Runge-Kutta  technique  is  a  multi-step  method  that  is  most  commonly  ap¬ 
plied  in  Computational  Fluid  Dynamics  or  Multi-Body  Dynamics.  As  an  explicit 
procedure,  it  is  computationally  efficient  and  the  myriad  of  schemes  provide  many 
alternatives  in  the  tradeoff  between  accuracy  and  speed.  Runge-Kutta  methods  are 
well-known  and  covered  extensively  from  numerous  viewpoints,  including  differential 
equations  [48,  142],  numerical  methods  [71,  77,  3,  103,  117],  and  computational  fluid 
dynamics  (CFD)  [14,  55].  Examples  of  its  application  with  conventional  FEA  are 
few  [94,  107].  However,  the  Runge-Kutta  method  is  commonly  employed  along  with 
the  Discontinuous  Galerkin  Method  [85,  53].  The  Runge-Kutta  method  discretizes 
the  first-order  differential  equation  in  time  to  generate  explicit  functions  for  slope  in 
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the  form  AU / At  =  H(U).  The  slope  is  calculated  at  optimized  points  within  the 
time  step.  Weighting  coefficients  refine  the  final  prediction  of  U  at  the  end  of  the 
time  step. 

Symplectic  Methods. 

The  symplectic  methods,  also  called  geometric  or  variational  integrators,  retain 
the  relationship  between  the  generalized  coordinates  q  and  their  generalized  velocities 
or  momenta  p.  Thus  the  flow  or  mapping  of  the  system  from  an  initial  point  (q,p) o  to 
its  future  configuration  (q,p)  1  should  preserve  the  area  of  a  given  set  of  points  as  they 
move  from  t  —  to  to  t  —  ti  [48].  In  a  discretized  Hamiltonian  system  the  coordinates 
and  momenta  are  related  to  the  Hamiltonian  as  Ap/At  =  —AH(q,p,t)/Aq  and 
Aq/At  =  AH(q,p,t)/Ap.  Often  in  mechanics  the  potential  energy  is  a  function 
only  of  position  ( U  =  U(q))  and  the  kinetic  energy  is  a  function  only  of  velocity 
( T  —  T(p)).  A  variational  integrator  iterates  between  q  and  p  across  multiple  stages 
to  achieve  a  prediction  for  the  future  configuration  vectors  q  and  p.  Because  the 
relationship  between  q  and  p  is  largely  preserved,  symplectic  integration  has  been 
shown  to  more  effectively  conserve  a  system’s  energy  over  time.  This  feature  is  of 
utmost  importance  during  long-duration  simulations  where  small  rates  of  energy  gain 
or  loss  result  in  substantial  errors  at  the  end. 

For  a  thorough  treatment  of  symplectic  integration,  refer  to  [48,  137,  84],  Also,  a 
seminal  work  in  the  held  of  symplectic  integration  is  that  of  Yoshida  [144] .  Yoshida 
developed  a  rigorous  method  for  deriving  exact  coefficients  for  high  order  integrators. 

Symplectic  integration  has  been  successfully  integrated  into  FE  models  of  beam 
vibration.  Leung  and  Mao  [82]  demonstrated  the  use  of  a  symplectic  algorithm  in 
the  FEA  of  non-linear  vibrations  of  a  beam  and  recommended  it  for  use  in  plates  and 
shells.  They  separated  the  kinetic  and  potential  energy  equations  to  utilize  the  mid- 
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point  symplectic  technique  [22],  which  they  called  the  “time-centered  Euler  scheme.” 
Dash  et  al.  [30]  used  Yoshida’s  8th  order  scheme  and  observed  good  results  but  marked 
instability  above  a  critical  At.  Note  that  some  of  the  other  methods,  such  as  the 
Newmark  and  certain  Runge-Kutta  formulations,  can  be  shown  to  be  symplectic  [68, 
84,  137].  Hughes  mathematically  proved  that  for  an  unforced  undamped  response  the 
Newmark  trapezoidal  rule  exactly  conserves  the  total  energy  of  a  system  [58]. 

The  fact  that  a  method  is  symplectic  does  not  necessarily  guarantee  increased 
performance  or  suitability  for  a  given  problem.  Energy  preservation  may  indeed  be 
counter  to  a  need  for  algorithmic  dissipation.  Thus,  in  this  research,  symplecity  will 
be  examined  to  gain  further  understanding  of  proposed  methods,  but  will  not  be  a 
basis  for  their  formulation  or  an  objective  to  be  pursued. 

Space-Time  Finite  Elements. 

Space-time  finite  element  formulations  have  been  seen  as  a  viable  alternative  to 
direct  integration  for  several  decades  [9,  17,  100].  In  this  method,  the  elements  span 
not  only  space  but  also  a  discrete  time  increment.  The  most  frequently  encountered 
space-time  hnite  element  technique  is  the  time-discontinuous  Galerkin  (frequently 
abbreviated  as  TDG  [62])  method.  In  the  discontinuous  Galerkin  method,  element 
shape  functions  produce  different  values  of  the  field  variable  at  their  interfaces;  flux 
schemes  are  then  employed  to  resolve  the  ambiguity.  High-order  shape  functions  using 
a  variety  of  basis  types  are  permitted.  For  example,  in  a  wave  problem,  basis  functions 
that  are  themselves  solutions  to  the  wave  equation  have  provided  great  accuracy  while 
minimizing  computational  cost  [101].  In  a  different  but  equally  effective  approach, 
Hamilton’s  Weak  Principle  has  been  applied  to  develop  simple  elements  with  constant 
and  linear  shape  functions.  The  integration  could  be  performed  by  inspection,  and 
led  to  a  system  of  algebraic  equations  for  solving  a  variety  of  nonlinear  dynamic 
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problems  [4,  56]. 

Implementation  of  space-time  FEA  in  all  of  the  aforementioned  studies  involved 
obtaining  the  solution  for  one  slab  of  elements  at  a  time,  where  a  slab  is  defined  as 
the  set  of  elements  in  a  given  time  increment  to  to  ti  —  to  +  At.  The  nodal  values  at 
f i  were  then  used  as  initial  conditions  for  the  following  time  increment.  In  contrast, 
the  method  employed  during  this  research  for  modeling  linear  membranes  solved  the 
entire  space-time  domain  simultaneously. 

Hermite  Time  Interpolation. 

An  alternative  time  integration  method  was  developed  for  this  study:  Hermite 
polynomial  interpolation.  A  unique  form  of  collocation,  this  scheme  uses  derivatives 
at  the  end  points  to  provide  the  necessary  number  of  equations  rather  than  additional 
interior  points  [81].  A  common  example  of  Hermite  polynomials  in  finite  element  anal¬ 
ysis  can  be  found  in  simple  beam  elements,  where  displacements  and  rotations  at  the 
ends  of  the  beam  define  the  cubic  polynomial  shape  functions  [107,  28].  Higher-order 
Hermite  interpolations  are  possible,  but  they  demand  additional  end  point  function 
derivatives,  interior  nodes,  or  some  other  means  of  providing  enough  constraints  to 
define  a  unique  solution. 

For  initial  value  problems  in  general,  a  wide  variety  of  time-marching  methods  are 
defined  by  the  selection  of  a  polynomial  interpolant  and  a  corresponding  numerical 
integration  scheme.  These  include  linear  multi-step  methods  and  single-step  colloca¬ 
tion  methods  like  the  implicit  Runge-Kutta  methods  [48,  3].  The  ubiquitous  implicit 
midpoint  rule  also  falls  under  this  category  [116].  These  methods  rely  on  internal 
collocation  points  for  their  definition.  In  contrast,  Hermite  polynomials  can  be  cast 
in  a  two-point  form,  so  the  interpolation  is  defined  by  the  function  value  and  its 
derivatives  at  the  end  points  only.  The  term  “prolongation”  has  been  used  to  refer  to 
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the  use  of  additional  derivatives  at  the  end  points  to  provide  the  necessary  number  of 
equations  rather  than  additional  interior  points  [81].  A  common  example  of  Hcrmite 
polynomials  in  finite  element  analysis  can  be  found  in  simple  beam  elements,  where 
cubic  shape  functions  are  uniquely  defined  by  the  displacements  and  rotations  at  the 
ends  of  the  beam  [107,  28].  Higher-order  interpolations  are  possible,  but  they  demand 
additional  end  point  function  derivatives.  One  can  also  obtain  a  solution  by  replacing 
one  of  the  derivatives  by  a  time  integral  constraint,  called  an  eliminant  [77]. 

Previous  results  using  Hermite  time  interpolation  have  shown  excellent  accuracy 
relative  to  the  central  difference  and  trapezoid  methods.  Leok  et  al.  [81]  used  quintic 
interpolation  polynomials  in  their  prolongation-collocation  method.  The  variational 
integrators  were  derived  by  applying  Euler- Maclaurin  quadrature  to  the  Lagrangian 
formulation.  It  was  mathematically  proven  that  a  quintic  polynomial  with  the  applied 
method  provided  at  least  fourth-order  convergence,  which  was  further  confirmed  by 
numerical  experiments.  The  experiments  also  showed  that  energy  was  conserved 
despite  structured,  bounded  fluctuations,  as  expected  for  a  variational  integrator. 
Stability  was  not  addressed. 

Quartic  Hermite  time  interpolation  was  applied  in  the  structural  dynamics  con¬ 
text  by  Razavi  et  al.  [109]  and  compared  to  members  of  the  Newmark  family.  The 
derivation  assumed  constant  mass,  damping,  and  stiffness  matrices,  so  it  applied 
only  to  linear  systems.  The  external  load  was  assumed  to  vary  linearly  within  the 
time  step.  The  unique  solution  was  defined  by  satisfaction  of  equilibrium  equation 
at  the  beginning  of  the  time  step,  at  the  end  of  the  time  step,  and  on  average  over 
the  time  step  (he.,  the  time  integral  of  the  residual  vanishes).  The  method  was 
labeled  a  weighted  residual  method  due  to  the  use  of  the  residual  equation  as  an 
eliminant.  Numerical  experiments  using  a  linear,  unforced  oscillator  corroborated 
the  fourth-order  convergence  rate  of  [81].  The  response  also  showed  no  dissipation 
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and  minimal  period  distortion.  Two  regions  of  instability  were  stated  in  approximate 
terms:  0.51  <  A t/T  <  0.55  and  A t/T  >  1.24,  where  the  dimensionless  sampling 
frequency  A  t/T  is  the  time  step  size  divided  by  the  signal’s  period  T. 

Summary 

In  this  chapter,  the  unique  characteristics  of  membranes  have  been  introduced, 
and  numerous  techniques  for  predicting  their  behavior  have  been  summarized.  The 
discussion  illustrates  why  certain  techniques  were  selected  for  the  current  research. 
For  example,  spring-mass  models  were  rejected  because  they  clearly  have  significant 
difficulty  representing  continuous  structures.  The  linear  membrane  model  enabled 
investigation  of  the  simultaneous  time-continuous  Galerkin  technique  because  it  re¬ 
quired  fewer  degrees  of  freedom  than  a  nonlinear  model,  a  key  consideration  because 
the  method  produces  an  extremely  large  linear  system  to  solve.  Likewise,  the  dis¬ 
continuous  Galerkin  method  was  rejected  because  it  increases  the  number  of  degrees 
of  freedom.  For  the  nonlinear  membrane  model,  the  point  collocation  spatial  dis¬ 
cretization  combined  with  the  group  FE  formulation  potentially  provided  simplicity, 
efficiency,  and  accuracy.  The  Hermite  time  interpolation  method  was  anticipated 
to  offer  higher  accuracy  with  acceptable  stability  limits.  Combining  point  colloca¬ 
tion  and  Hermite  interpolation  into  a  dynamic  membrane  model  results  in  a  purely 
collocation-based  scheme,  a  unique  and  promising  formulation  for  a  membrane  prob¬ 
lem. 
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III.  Linear  Analysis  of  Membrane  Dynamics 


This  chapter  details  the  development  of  a  space-time  finite  element  (FE)  model 
to  simulate  a  linear  membrane.  Space-time  FE  formulations  of  non-periodic  systems 
typically  utilize  a  time-marching  algorithm.  The  solution  is  obtained  for  one  slab  of 
elements,  where  a  slab  is  defined  as  the  set  of  elements  in  a  given  time  increment  to 
to  ti  —  to  +  At.  The  nodal  values  at  1 1  are  then  used  as  initial  conditions  for  the 
following  time  increment. 

The  goal  of  this  portion  of  the  research  was  to  investigate  a  space-time  FE  method 
simultaneously  discretized  and  solved  across  the  entire  space-time  domain.  This  ap¬ 
proach  results  in  extremely  large  linear  systems.  Because  the  discontinuous  Galerkin 
methods  with  low-order  interpolations  require  more  degrees  of  freedom  for  the  same 
mesh  [145],  continuous  interpolation  functions  in  both  space  and  time  (trilinear  el¬ 
ements)  are  employed.  Thus,  using  the  terminology  of  [62],  it  is  labeled  a  time- 
continuous  Galerkin  (TCG)  method.  Further,  the  proposed  approach  will  be  referred 
to  as  the  simultaneous  time-continuous  Galerkin  (STCG)  method  to  distinguish  it 
from  time-marching  methods. 

After  deriving  the  governing  equations  for  a  mixed  formulation,  stability  will  be 
examined  through  numerical  experimentation,  and  static  and  dynamic  verification 
will  be  used  to  assess  the  accuracy  of  the  scheme. 

Membrane  Model 

The  linear  membrane  model  described  in  Section  II  is  employed,  which  classically 
leads  to  the  linear,  second-order  PDE  of  Eq.  1.  With  displacement  w(x,y,t),  specific 
density  p(x,y),  and  tension  per  unit  length  P(x,y),  the  kinetic  energy  density  (the 
kinetic  energy  per  unit  area,  or  over  the  infinitesimal  spatial  area  d.tt  where  hi  is  the 


spatial  domain)  is 


Tdo,  — 


(17) 


The  strain  energy  density  is  a  function  of  the  membrane  gradients  and  the  tension. 


Udn  — 


(18) 


Mixed  Form. 


Discretization  of  the  governing  equation  in  its  current  form  results  in  a  displace¬ 
ment  formulation.  To  obtain  a  mixed  formulation,  the  second-order  equation  will  be 
resolved  into  a  system  of  first-order  equations  by  the  application  of  constraints.  Using 
the  variational  approach,  the  governing  equations  will  now  be  re-derived  using  La¬ 
grange  Multipliers  to  introduce  new  dependent  variables  and  enforce  the  constraints. 

The  generalized  momentum  p  and  generalized  stresses  a  are  defined  as 


P  = 

= 


dTdn  _  _dw 
dw  P  dt 
dUdn  =  pdw 
dw/dx  dx 
dUdn  _  pdw 
dw/dy  dy 


(19) 

(20) 

(21) 


The  distinction  that  these  terms  are  “generalized”  indicates  the  variables  are  not 
necessarily  based  upon  a  physical  model.  In  particular,  the  generalized  stresses  should 
not  be  confused  with  Cauchy  stresses.  These  relationships  will  be  enforced  by  the 
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Figure  1.  Diagram  of  degrees  of  freedom  for  the  mixed  formulation. 


constraint  equations  g^  =  0,  where 


9i=P~  pw 


92  =  ~  P 

93  =  oy-  P 


dw 

dx 

dw 

dy 


(22) 

(23) 

(24) 


Taking  the  variation  of  the  constraint  equations  for  use  later  in  the  derivation,  we 
obtain 


Sgi  = 
5g2  = 
5 93  = 


6p-  p 


d5w 


4  (T) 


dt 

IDdSw 

dx 
ddw 

-  p^r 


(25) 

(26) 
(27) 


In  terms  of  the  generalized  momentum,  the  variation  of  the  kinetic  energy  expression 
becomes 

h7do  =  piuSw  =  -p5p  (28) 

P 
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Similarly,  in  terms  of  generalized  strains,  the  potential  energy  becomes 


Udn  =  (a2x  +  ciy)  (29) 

and  its  variation  is 

SUdd  =  ~p  (°xScrx  +  <Jy5ay )  (30) 

The  virtual  work  due  to  a  prescribed  external  pressure  force  /,  where  in  the  linear 
membrane  case  f(x,  y,  t )  acts  in  the  same  direction  as  the  displacement  w,  is  expressed 
as 

=  fSw  (31) 

Applying  Hamilton’s  Law  of  Varying  Action. 

To  apply  Hamilton’s  Law  of  Varying  Action,  the  energy  functional  is  formed  . 

/  -  J  i^T  -  U  +  \tg)j  dt  (32) 

Next,  the  sum  of  virtual  work  and  the  first  variation  of  the  energy  functional  is  set 
to  zero  [7]. 

51  +  5W  =  0  (33) 

The  expressions  for  the  individual  terms  are  substituted,  and  the  chain  rule  is  applied 
to  the  Lagrange  Multiplier  summations,  leading  to  the  expanded  expression 

t\  f 

0  =  |  JJ  -5p  -  pS(7x  -  °pday  +  ^  A  iSg.i  +  ^  grfXi  +  fSw  dttdt  (34) 

t0  n  L  ii 

Substitution  of  the  constraint  functions  results  in  the  functional 
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-Sp  -  °^8ax  -  ^ 5ay  +  f5w 
p  P  P 


+  Ai  [Sp-p- 


+  Ao  (  5ax  —  P- 


+  A3  [  5(7 v  —  P- 


+  (p  -  P^j  6Xi  +  (^x  -  P^j  SX2  +  (ay  -  P^Pj  6X3  dPldt  (35) 


Integration  by  parts  is  applied  to  each  term  that  contains  a  derivative  of  a  variation 
of  a  variable.  This  procedure  moves  the  derivatives  to  the  Lagrange  multipliers  and 
extracts  the  boundary  conditions. 


/  /  /  p\1^——dUdt=  I  I  f  p^-P-5w  dVLdt  —  I  f  pXiSwdPL 


to 


to  Q 


PXo—dttdt  = 


I II  pd^6wdndt- 1 J 

to  ft  to  r 


PX2SwnxdTdt  (37) 


ti  t\  t 1 

j  jj  PX3-^—dYldt  —  I  J  J  P-^-8wy  dYldt  —  J  j  PX38wnydTdt  (38) 

t0  n  t0  n  to  r 
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After  performing  the  integration  by  parts,  the  functional  is 


0  = 


a. 


(7  ii 


-dp  ~  ~^dax  ~  ~^day  +  f$w 
p  P  P 


to  fl 


.  _  .  _  .  _<9Ai .  <9 A 2  f  _9A 3  r 

+  Ai  dp  +  A  28<jx  +  \3<jy  +  p^—6w  +  P——5w  +  P— — <5w 

at  ox  oy 

_dw\  (  dw A  /  ndw\ 

P-PW)  «1+U-P^ ) 


dildt 


p\i  5w  dPt 


Q 
1 1 


[P A2  die  nx  +  P\3dw  ny\  dr dt 


(39) 


to  r 


Next  we  collect  terms  associated  with  each  variational  term. 


t0  r 
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Since  this  equation  must  be  true  for  all  permissible  functions  Sw ,  5p,  5ax,  and  Sayi 
the  meanings  of  the  Lagrange  multipliers  can  be  determined  by  inspection. 


Ai 

A2 

A3 


_P_ 

P 

O X 

~p 

ay 

P 


(41) 

(42) 

(43) 


Substituting  these  expressions  into  the  functional  eliminates  the  Lagrange  multipliers. 


to  r 

By  collecting  the  terms  according  to  their  variational  components,  we  obtain  the 
mixed  form  of  the  governing  equation  and  recover  the  three  constraint  equations. 
These  four  equations  collectively  form  the  basis  for  the  upcoming  finite  element  de¬ 
velopment. 


dp 

dt 


do x  _  day 
dx  dy 
_dw 


P~  P 
<TX—P 

Gy~P 


dt 

dw 

dx 

dw 

dy 


f 

0 

0 


(45) 
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0 


(46) 

(47) 

(48) 


Domain  Discretization 


The  space-time  mesh  was  constructed  by  duplicating  the  spatial  mesh  at  each 
time  increment  t  +  nAt.  Each  initial  spatial  node  (x,y,to),  referred  to  in  this  paper 
as  a  parent  node,  has  the  offspring  ( x,y,t0  +  nAt).  Thus  each  node  in  the  domain  is 
uniquely  identified  by  its  ( x ,  y,  t)  coordinates.  For  this  study  At  was  constant,  though 
it  could  certainly  be  varied  to  optimize  performance.  Each  node  has  four  degrees  of 
freedom:  w,  p,  crx,  and  ay,  any  of  which  could  be  fixed  or  otherwise  prescribed. 

Elements. 

Two  well-known  three-dimensional  solid  elements  were  utilized  as  space-time  el¬ 
ements.  Both  are  basic  trilinear  isoparametric  elements.  Using  £  and  r/  to  denote 
in-plane  local  coordinates,  the  local  coordinate  system  is  denoted  as  (£,  r/,  r)  and  is 
aligned  with  the  global  system  (. x,y,t ).  Though  the  third  axis  is  time  rather  than 
space,  no  modifications  to  the  interpolation  functions  or  integration  procedures  were 
necessary.  Conceptually  the  resulting  elements  may  be  thought  of  as  two  planar  mem¬ 
brane  elements,  one  at  time  t  and  the  other  at  time  t  +  At,  stacked  upon  each  other 
to  create  an  internal  volume  that  is  a  space-time  subdomain  T>.  A  slab  is  defined 
as  the  set  of  all  of  the  elements  in  space  between  two  given  times  [60],  so  the  entire 
solution  domain  consists  of  a  stack  of  Nt  slabs  where  Nt  is  the  number  of  time  steps. 
In  this  paper,  the  term  “surface  element”  refers  to  a  quadrilateral  or  triangle  face  in 
the  (x,  y)  plane  that  is  shared  between  two  elements  at  the  given  time  slice. 

The  first  element  is  variously  called  a  trilinear  hexahedral  [58],  hexahedron  [102], 
or  when  using  right  angles  a  linear  brick  [58].  It  is  an  extension  of  the  bilinear 
quadrilateral  into  the  third  dimension.  The  complete  interpolation  polynomial  is 
Nt  —  (1  +  £;£)(1  +  ryy)(l  +  7iT)/ 8  where  £,y,  r  e  [—1,1],  providing  C°  continuity 
between  elements.  The  4th-order  accurate  2x2x2  integration  points  are  (£,  rj,  r)  = 
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(a)  Hexahedra  (b)  Triangular  prisms 

Figure  2.  Representative  slabs  of  elements. 

(Ay/l/3),  (±y/l/3),  (Ay/l/3).  Reduced  integration  performed  poorly  and  was  not 
used. 

The  second  element  is  the  triangular  prism,  also  called  a  wedge  or  a  right  penta¬ 
hedron  [102].  It  can  be  viewed  as  an  extension  of  the  constant-strain  triangle  or  a 
degenerate  case  of  the  hexahedron  [58].  The  C°  linear  shape  functions  are  given  in 
Eq.  49  [67]  with  £,77  e  [0, 1]  and  r  G  [—1, 1]. 


Ni  =  ^(l-£-r])(l  +  TiT),  i  =  l,4  (49) 

Ni  =  ^£(1  +  rt),  i  =  2,5 
Ni  =  ^77(1  +  rt),  7  =  3,6 

The  3x2  integration  points  are  (£,77)  =  (|,  |),  (|,  |),  (|,  |)  at  r  =  ±y/ 1/3,  suffi¬ 
cient  for  up  to  a  second  degree  polynomial.  As  with  the  hexahedral  element,  reduction 
of  the  integration  order  severely  degraded  the  results. 

The  local  element  matrices  were  constructed  to  satisfy  the  governing  equations  as 
shown  in  Eq.  50.  Since  there  are  four  degrees  of  freedom  per  node  the  resulting  block 
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stiffness  matrix  K  was  32x32  for  a  hexahedron  and  24x24  for  a  triangular  prism. 
Boundary  conditions  were  imposed  locally  by  removing  the  appropriate  rows  and 
columns  from  the  stiffness  matrix  and  adjusting  the  right-hand  side  accordingly  [28]. 
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(51) 

(52) 

(53) 

(54) 


where  T>  is  the  element  subdomain  of  the  problem  domain  12  x  [to,tf\.  The  non¬ 
zero  values  were  then  inserted  into  the  sparse,  square  global  stiffness  matrix.  For 
a  discretization  with  Ns  spatial  nodes  and  Nt  time  steps,  and  accounting  for  the 
four  degrees  of  freedom  per  node  due  to  the  mixed  formulation,  the  matrix  size  is 
[4 Ns(Nt  +  l)]2  minus  rows  and  columns  eliminated  through  imposition  of  boundary 
or  initial  conditions.  For  dense  spatial  meshes  or  long  simulations,  hardware  ca¬ 
pacity  clearly  becomes  a  practical  concern  and  at  some  point  this  method  becomes 
unsuitable.  The  formulation  results  in  a  single  linear  equation,  the  solution  of  which 
provides  the  displacement,  momentum,  and  in-place  stresses  for  each  node. 

After  obtaining  the  solution  vector,  post-processing  matched  the  degrees  of  freedom  to 
the  parent  nodes  to  create  the  response  history  for  each  surface  node  of  the  membrane. 
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Static  Verification  and  Spatial  Convergence. 


For  the  present  method,  a  static  problem  is  a  special  case  of  the  general  method  - 
the  same  algorithm  is  used  for  both  static  and  dynamic  analyses.  For  a  static  analysis 
the  generalized  momentum  p  is  set  to  zero  for  every  node  in  the  space-time  mesh,  and 
spatial  boundary  conditions  are  fixed  to  their  values  at  t  —  0.  These  constraints  are 
sufficient  to  determine  a  unique  static  solution.  Though  only  one  time  step  is  required 
to  generate  a  single  slab  of  elements  and  subsequent  static  solution,  the  solution  is 
identical  for  an  arbitrarily  large  number  of  time  steps.  The  equilibrium  of  each  static 
solution  was  verified  during  post-processing  by  checking  that  the  field  values  of  a 
parent  node  were  identical  to  those  of  its  offspring. 

Static  verification  was  performed  by  use  of  Examples  2.8  and  4.1  in  Reddy’s  text¬ 
book  [110],  which  conveniently  provide  both  an  analytical  solution  and  finite  element 
approximations.  The  FE  solutions  were  obtained  using  displacement  formulation  and 
bilinear  quadrilateral  and  triangle  elements  on  two  grids  (element  size  he=  0.25  and 
he  =  0.125,  or  equivalently  as  the  number  of  elements  in  his  mesh,  N  =  4  and  N  =  8). 
Setting  P  —  p  —  f  —  1,  the  governing  equation  is  the  Poisson  equation, 

-  V2w  =  1  (55) 

The  domain  was  a  1  x  1  square,  x,y  e  [0, 1],  with  mixed  boundary  conditions 


w(x,  1)  =  0 

(56) 

w(l,y)  =  0 

(57) 

8m 

or(°^)=° 

(58) 

8m 

(59) 
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The  series  solution  was  therefore 


(60) 


The  STCG  method  was  employed  using  both  types  of  elements.  The  error,  e,  was 
defined  as  the  displacement  error  at  the  origin.  The  observed  rate  of  convergence  for 
both  elements  was  two,  as  shown  in  Figure  4.  The  errors  were  approximately  the 
same  as  those  of  Reddy’s  displacement  method. 
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Figure  3.  Representative  Poisson  equation  solution  using  square  elements  ( he  =  0.1). 

Figure  5  plots  STCG  results  on  the  x-axis  for  the  two  coarsest  mesh  sizes  (the  same 
as  used  by  Reddy).  Mild  variations  of  the  solution  are  evident,  a  phenomenon  that 
can  occur  in  mixed  formulations  with  certain  interpolation  combinations  [19,  146].  In 
this  case  they  reduced  significantly  as  the  mesh  was  refined  and  the  model  remained 
consistent. 

Time  Integration  Performance. 

In  this  section  we  investigate  the  accuracy  and  reliability  of  the  model  in  repro¬ 
ducing  periodic  signals  of  varying  frequency.  The  objective  is  to  derive  a  methodology 
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for  determining  At  and  characterizing  limitations  and  restrictions  on  the  choice.  At 
least  ten  samples  per  period  are  required  by  rule  of  thumb  [120]  and  as  the  practical 
limit  for  a  linear  approximation  of  a  periodic  signal  [130].  These  guidelines  as  well  as 
stability  limits  will  be  tested  by  numerical  experimentation  and  compared  to  other 
methods.  Let  to  be  the  angular  frequency.  Frequencies  will  be  normalized  by  the 
sampling  rate  oos.  Depending  on  context,  the  normalized  frequency  will  be  written 
either  as  c o*  where  to*  =  oo/ooS)  or  as  A t/T,  which  is  numerically  identical  and  found 
in  finite  element  texts  [8,  58].  To  avoid  high  frequencies  being  folded  into  the  lower 
frequency  band  (aliasing)  one  must  sample  at  more  than  twice  the  Nyquist  frequency 
(cu;  >  2,to n,  where  t o%-  =  0.5)  [105,  33,  120].  Finally,  since  signals  are  distorted  by 
the  model,  frequencies  observed  in  model  output  are  denoted  with  a  bar  (e.g.,  a)*).  It 
should  be  noted  that  a  spectrum  of  normalized  frequencies  is  present  in  every  dynamic 
response  signal. 

First  we  look  at  how  the  sampling  rate  affects  the  Gaussian  quadrature  accuracy  in 
the  time  axis.  Consider  the  time  span  of  one  element  slab  At  and  a  periodic  signal  of 
period  T  traversing  it.  The  fraction  of  the  period  covered  by  the  element  has  already 
been  defined  as  the  normalized  frequency  to*  =  A  t/T.  Consider  also  that  the  element 
may  be  sampling  any  portion  of  the  signal  depending  on  the  phase  <fi.  The  error  is 
the  percentage  difference  between  the  results  of  the  well-known  two-point  Gaussian 
quadrature  integration  and  the  analytical  solution.  The  bounds  of  the  integration 
error  shown  in  Figure  6  depict  the  normalized  error  at  the  worst-case  (j).  As  the 
sampling  rate  approaches  infinity  (A t/T  — >  0),  the  error  approaches  zero  because 
the  linear  approximation  approaches  the  exact  solution  regardless  of  the  phase.  Also 
regardless  of  phase,  when  sampling  at  the  Nyquist  frequency  the  linear  quadrature 
produces  a  zero  result  (the  integration  points  will  be  equidistant  from  the  t-axis  and 
on  opposite  sides).  The  error  is  less  than  3%  for  to*  =  0.1  and  less  than  13%  for 
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uj*  =  0.2. 


Figure  6.  Worst-case  quadrature  error  bounds  for  all  phase  angles. 


For  the  following  numerical  investigations  a  single  degree-of-freedom  simple  har¬ 
monic  oscillator  was  formed  using  two  equally-sized  surface  elements  in  the  (x,  y) 
plane.  Their  edges  were  joined  along  the  axis  x  —  0,  they  were  placed  between  simple 
supports  located  at  x  =  ±1,  subjected  to  an  initial  displacement  tc(0,?/,0)  =  1  along 
their  interface,  and  released.  The  solution  was  w  =  cost ct  where  u  =  1.5c/ h. 


Initial  displacement  w  =  1, 
then  released 


Figure  7.  Diagram  of  the  simple  harmonic  oscillator  configuration. 


First  a  convergence  study  in  time  was  performed.  The  simulation  was  run  for 
two  periods  with  varying  numbers  of  elements  (equivalent  to  varying  the  normalized 
frequency).  The  error  was  defined  as  the  mean  error  magnitude  for  all  solution  points. 
As  shown  in  Figure  8,  the  observed  rate  of  convergence  was  two. 

Next  the  oscillator  was  run  for  100  periods  over  a  series  of  runs  with  different  values 
of  oj*  for  the  oscillator.  The  values  of  u f  were  varied  by  changing  the  oscillator’s 
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Figure  8.  Time  convergence  results  for  the  STCG  model. 

natural  frequency  with  E.  As  seen  in  Figure  9,  a  spurious  mode  was  evident  at 
a  higher  frequency  than  the  fundamental  response.  As  u*  moves  right  from  the 
origin  the  spurious  mode  migrates  leftward  from  ujn-  This  behavior  often  indicates 
the  spurious  peak  has  been  aliased.  The  frequencies  of  the  peaks  given  the  known 
input  frequency  are  shown  in  Figure  10(a).  Given  an  to*  on  the  vertical  axis,  one 
reads  horizontally  right  to  find  the  model-predicted  normalized  frequencies  uj*  of 
both  responses.  In  every  case  the  peaks  merged  when  c o*  ~  0.276  and  uj*  ~  1/3.  The 
amplitude  distortion  as  the  frequency  increases  may  be  seen  in  Figure  10(b)  where 
the  largest  observed  displacement  and  the  mean  of  the  displacement  magnitudes  for 
the  entire  run  are  depicted  for  various  runs.  For  a  perfectly  reproduced  signal,  the 
maximum  would  be  1.0  and  the  average  magnitude  approximately  0.64.  A  ramp-up  of 
the  maximum  observed  displacement  was  evident  prior  to  the  peaks  merging.  When 
the  frequency  of  the  oscillator  was  raised  further  the  solution  was  severely  damped 
and  converged  to  zero  within  a  few  times  steps.  Though  the  precise  mechanism  in  this 
model  has  not  yet  been  determined,  this  behavior  has  been  noted  in  finite  element 
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complex  wavenumber  Fourier  analyses  [130].  Accordingly  we  will  refer  to  a)*  ~  0.333 
as  the  cutoff  frequency. 

Nearly  identical  behavior  was  mathematically  detailed  and  reproduced  by  [130] 
while  performing  a  finite  element  complex  wavenumber  Fourier  analysis.  The  fol¬ 
lowing  briefly  describes  their  findings  for  comparison  to  the  present  results.  Linear 
element  discretization  of  the  reduced  wave  equation  for  a  uniform  elastic  bar  produces 
regions  of  the  frequency  domain  called  the  passing  band  and  the  stopping  band.  They 
are  separated  by  the  cutoff  frequency  amax  using  a  different  normalization  a  =  uh/c. 
When  a  <  amax  the  imaginary  component  of  the  eigenvalue  is  zero  and  the  wave 
propagates  (passes);  when  a  >  amax  the  imaginary  components  produce  rapid  damp¬ 
ing,  as  seen  in  our  Figure  10(b).  For  a  diagonalized  mass  matrix  amax  =  2,  and 
the  reference  frequency  aref  is  set  at  ten  elements  per  wavelength  (a/  =  0.1  in  our 
nomenclature).  Then  we  have  from  their  development  the  ratio 


Q-max 

Q-ref 


2 

2tt/10 


10/ 7T 


(61) 


Setting  the  ratio  h/c  constant,  the  ratio  amax/aref  is  equal  to  the  ratio  of  equivalent 
normalized  frequencies  in  our  development.  Hence  we  can  solve  for  our  cutoff  fre¬ 
quency.  For  the  simple  harmonic  oscillator  uj*  —  0.1  and  c os  =  15  rad/s  produced  an 
output  signal  at  to*  =  1.52  rad/s. 


—  * 


0.323 


(62) 


This  prediction  is  within  3%  of  our  model’s  cutoff  frequency.  The  comparison  between 
results  strongly  suggests  that  the  current  technique  is  filtering  out  all  content  above 
a;*  «  0.276. 

The  numerical  dispersion  in  terms  of  period  distortion  was  compared  to  that  pro- 
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(a)  Displacement,  w*  =0.1 


(c)  Displacement,  w*  =  0.2 


Normalized  frequency 

(b)  Frequency,  w*  =0.1 

. 

100 


Normalized  frequency 

(d)  Frequency,  w*  =  0.2 


(e)  Displacement,  w*  =  0.27  (f)  Frequency,  w*  =  0.27 


(g)  Displacement,  w*  =  0.28  (h)  Frequency,  oj*  =  0.28 


Figure  9.  Displacement  history  and  associated  frequency  response  of  the  simple  har¬ 
monic  oscillator  with  increasing  time  steps. 
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(a)  Frequencies  of  true  and  spurious  modes. 
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(b)  Amplitude  distortion  with  increasing  u;* . 


Figure  10.  Response  frequency  and  magnitude  behavior  of  the  simple  harmonic  oscil¬ 
lator. 


duced  by  two  of  the  most  commonly-used  direct  integration  techniques:  the  central 
difference  method  and  the  trapezoidal  rule.  Both  are  sub-cases  of  the  Newmark 
method.  The  central  difference  method  is  obtained  when  (3  =  0  and  7  =  1.  It  is 
explicit  and  therefore  conditionally  stable  with  the  limit  A t/T  <  1/n.  The  trape¬ 
zoidal  rule,  where  (3  =  \  and  7  =  \ ,  is  symplectic  with  unconditional  stability.  The 
dispersion  relationship  for  the  present  method  can  be  seen  on  the  left-hand  curve  of 
Figure  10(a),  which  shows  a  nearly  linear  relationship  of  oj*/Cj*  through  u*  =  0.2. 
Separately,  the  one-dimensional  central  difference  and  trapezoid  rule  algorithms  were 
used  to  solve  the  same  ordinary  differential  equation  [102],  The  results  are  shown 
in  Figure  11.  As  the  normalized  frequency  increased,  the  trapezoid  rule  elongated 
the  period  and  the  central  difference  method  compressed  it,  matching  results  in  the 
literature  [8].  The  simultaneous  TCG  preserved  the  period  to  within  0.1%  up  to 
A  t/T  =  0.1  and  within  0.7%  up  to  A  t/T  =  0.18.  This  frequency  preservation  is  vital 
since  our  primary  goals  include  the  reproduction  of  modal  natural  frequencies. 

Although  the  present  model  was  not  susceptible  to  instability,  the  cutoff  frequency 
places  a  clear  limit  on  the  useful  part  of  the  spectrum.  From  the  cutoff  frequency 
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Figure  11.  Period  distortion  comparison  for  different  integration  schemes 

one  can  derive  the  maximum  allowable  time  step  A tcr  in  terms  of  element  size  and 
wave  speed.  Note  that  the  standard  notation  h  is  used  here  for  the  representative 
maximum  element  spatial  length,  not  to  be  confused  with  membrane  thickness. 


(63) 


The  constant  1.16  represents  the  maximum  allowable  Courant  Friedrichs-Lewy  Num¬ 
ber  CFL.  The  cutoff  frequency  allows  a  slightly  higher  CFL  than  the  broadly  rep¬ 
resentative  cases  listed  in  Table  1  and  is  therefore  less  restrictive  on  the  choice  of 
At. 

Table  1.  Survey  of  representative  values  of  CFLmax  from  the  literature. 


1.0 

Wave  equation  with  linear  elements  [58] 

1.0 

Unsupported  axial  bar  using  lumped  masses  [28] 

0.707 

Rectangular  4-node  quadrilateral  elements  with  central  difference  method  [58] 

0.577 

Two- node  linear  rod  with  central  difference  method  [58] 

When  solving  the  wave  equation  with  a  linear  element  and  the  Newmark  method, 
CFL  =  1  is  not  only  limiting  but  also  optimal  [58].  Similarly,  for  a  two-node  linear 


47 


rod  element  using  the  Newmark  central  difference  method,  the  limit  is  CFL  =  0.577. 


At  <  «  0.577-  (64) 

V3c  c 

Hughes  also  provided  the  results  of  studies  of  rectangular  4-node  quadrilateral  ele¬ 
ments  along  with  the  central  difference  method,  which  resulted  in  CFL  =  0.707: 


At  < 


cd(l/h\  +  1  /hi) 


1/2 


h 

0.707- 

c 


(65) 


where  for  generality  we  equalize  the  spatial  dimensions  h  «  hi  «  h2  and  recall  that 
cj  =  (\  +  2/i)/ p  =  Ej p  as  found  earlier.  Lastly,  Cook  [28]  modeled  the  mechanism 
as  an  unsupported  axial  bar  of  length  h.  The  highest  frequency  is  umax  =  irc/h  since 
masses  are  not  lumped.  Thus  he  obtains  CFL  =  0.637  and  the  critical  time  step 


A tcr  < 


2 


^ max 


2  h 

TIC 


h 

0.637- 

c 


(66) 


Conservation  of  energy  was  examined  by  running  the  simple  harmonic  oscillator 
for  100  periods.  This  time  scale  was  sufficient  for  this  study  since  the  dynamic  cases 
capture  at  most  five  periods  of  the  first  mode.  The  phase  plots  for  oscillations  at 
three  different  normalized  frequencies  are  shown  in  Figure  12.  At  co*  =  0.1  the  energy 
was  effectively  conserved.  At  u*  =  0.2,  where  integration  error  up  to  13%  has  been 
demonstrated,  the  plot’s  slow  migration  towards  the  origin  indicates  a  progressive  loss 
of  energy.  At  oj*  =  0.28,  just  above  the  cutoff  frequency,  the  extreme  energy  dissipa¬ 
tion  was  obvious  as  the  solution  collapsed  to  the  resting  position.  This  phenomenon 
was  also  clear  in  terms  of  amplitude  decay  in  Figure  10(b). 

In  this  section  we  have  established  the  following  performance  observations:  fre¬ 
quencies  are  accurately  reproduced  through  c o*  =  0.10,  accuracy  degrades  through 


48 


Figure  12.  Energy  conservation  of  the  simple  harmonic  oscillator  at  different  normalized 
natural  frequencies. 


c o*  =  0.20  clue  to  integration  error,  and  results  are  increasingly  unreliable  as  the 
cntoff  frequency  uj*  ~  0.276  is  approached.  Since  the  model  aggressively  dissipates 
signals  above  the  cntoff  frequency,  frequency  content  there  is  noise.  The  next  section 
outlines  a  methodology  based  on  these  results. 


Membrane  Dynamics 

Cases  and  Methodology. 

Two  canonical  cases  were  selected  for  comparing  model  results  to  analytical  solu¬ 
tions:  a  rectangular  membrane  and  a  circular  membrane.  The  setup  parameters  are 
provided  in  Table  2.  Note  that  only  one-fourth  of  the  membrane  was  modeled  in  each 
case  using  zero-force  boundary  conditions.  The  steps  of  the  method  follow: 


1.  Choose  the  target  frequency  cutgt ,  the  highest  frequency  that  must  be  accurately 
reproduced. 

2.  Calculate  the  time  increment.  To  ensure  sufficient  accuracy  at  cotgt, 


At 


Ttgt 

10 


71 

5u>tgt 


(67) 
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3.  Ensure  the  surface  elements  are  properly  sized  (h  is  element  size).  From  Eq.  63, 


h  >  - =  0.865cAtcr  (68) 

O  r  Lsynax 

As  noted  in  the  literature  for  wave  equation  discretization  [58,  28],  CFL  values 
close  to  1.0  provided  the  best  results.  When  CFL  was  near  0.5,  for  example, 
oscillations  were  introduced  at  wavelengths  that  interfered  with  the  targeted 
modes.  Thus,  for  best  results  strive  for 

h  =  cAt  (69) 


4.  Run  the  simulation.  All  nodes  at  t  —  0  require  initial  displacements  and  mo¬ 
menta. 


5.  Determine  the  highest  valid  frequency  content  in  the  model’s  results. 

(a)  As  demonstrated  by  the  cutoff  frequency,  content  with  T  <  3A t  is  not 
valid  (the  small  wavelengths  are  not  sampled  at  a  sufficient  rate). 

(b)  When  the  CFL  is  low  (element  size  large  relative  to  the  time  increment), 
the  oscillations  span  more  time  intervals.  From  Eq.  63  and  the  size  limit 
in  Eq.  68,  the  noise  region  is  defined  as  T  <  (47t/3 ){h/c). 

After  the  results  of  the  model  have  been  obtained,  one  may  consider  invalid  any 
content  with  a  period  of 


T  <  max 


(70) 


The  analytical  solution  for  the  rectangular  membrane  is  given  in  Eq.  71  [108]. 
To  target  the  third  modal  frequency,  At  was  calculated  to  be  one  tenth  of  the  third 
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mode’s  wavelength.  The  membrane  was  subjected  to  an  initial  velocity  that  was 
uniform  across  the  entire  surface. 

Table  2.  Dynamic  test  cases  for  the  STCG  model.  Only  one  quadrant  was  modeled 
(x  >  0  and  y  >  0). 


Rectangle 

Circle 

Forcing  function 

Uniform 

r  <  0.3a:  Impulse  p0  in  [t0,ti] 

initial  velocity 

r  >  0.3a:  Zero 

Element  type 

Right  Hexahedra 

Triangular  prisms 

Dimensions 

2.4  x  2.0 

unit  radius(a  =  1) 

Surface  mesh 

143  nodes 

81  nodes 

120  elements 

121  elements 

Element  length  h 

0.141 

0.110 

c 

1.0 

1.0 

At 

0.128 

0.114 

CFL 

0.91 

1.04 

Stiffness  matrix 

48,454  square 

31, 424  square 

2,  305,  821  non-zero 

1, 193,  548  non-zero 

0.10%  non-zero 

0.12%  non- zero 

wmn{x,y,t)  = 


mirx  .  niry 


sm 


sm 


m= 1  n= 1 


a 


-R \nn  sill  jJrnn t 


(71) 


^ m.n  —  C7T 


Bmn 


/ra\2  / n 

Vo/  V  b 


a  nb 


1/2 


abu . 


mirx  nny 

w o(x,  y)  sm - sm  — - — axay 

run  JO  JO  ®  0 

4 

(cOSm7T  —  1)  (cOS  T17T  —  1) 


mrni^Lumn 


For  the  circular  membrane  the  two-dimensional  mesh  as  shown  in  Figure  13(b)  was 
generated  using  the  open  source  software  Gmsh  [43].  The  membrane  had  an  outer 
radius  a  =  1  and  the  impulse  was  a  uniform  pressure  p0  applied  inside  the  radius 
b  =  0.3  only  at  t  =  0.  The  natural  frequencies  were  umn  =  c/3mn  where  f3mna  are  the 
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zeros  of  the  Bessel  functions  Jm(f3a).  The  exact  solution  is  given  in  Eq.  72  [91].  The 
time  increment  was  chosen  to  target  the  third  modal  frequency. 


(a)  Rectangular  membrane.  (b)  Circular  membrane. 

Figure  13.  Meshes  used  for  evaluation  of  the  STCG  method  for  the  dynamic  cases. 


w(r ,  9,  t ) 


2p06c  y-v  Jl  [{UQn/c)  b\  J0  [(cjp n/c)  T 
pa2  [(^On/c)  a] 


f{r)sinu 0n  (t 


t)  dr 


(72) 


Modal  Frequency  Results. 

The  power  spectral  density  (PSD)  for  the  rectangular  membrane  with  respect  to 
co*  for  every  node  is  shown  in  Figure  14(a).  The  vertical  lines  indicate  the  exact 
solutions  co*  so  any  phase  distortion  can  be  discerned  by  offset  peaks.  The  first  seven 
mode  frequencies  were  accurately  predicted  but  beyond  do  =  8  rad/s  (approximate 
co*  =  0.16)  the  calculated  modal  frequencies  are  inaccurate.  Figure  14(b)  displays  a 
single  PSD  that  is  the  supremum  of  all  nodal  values  at  each  frequency;  thus  all  PSDs 
in  Figure  14(a)  lie  under  the  single  PSD  in  Figure  14(b).  The  dotted  line  indicates  the 
exact  solution  for  the  discrete  nodes  and  was  thus  equally  affected  by  the  sampling 
rate.  The  first  seven  modes  were  included  in  the  exact  solution. 

The  PSD  plots  for  the  circular  membrane  were  created  in  the  same  manner  as  those 
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10° 


Angular  Frequency  (rad/s) 


(a)  Composite  plot  of  PSDs  for  every  node.  (b)  Comparison  of  model  results  to  exact  so- 

Vertical  lines  depict  exact  modal  frequencies.  lution. 


Figure  14.  Frequency  spectrum  of  the  response  of  the  rectangular  linear  membrane 


of  the  rectangular  membrane  and  are  shown  in  Figure  15.  The  analytical  solution  in 
Figure  15(b)  includes  only  the  first  four  modes.  The  first  three  mode  frequencies  were 
accurately  reproduced  with  the  third  mode  at  approximately  to*  =  0.15.  However, 
the  spatially  discontinuous  impulse  generated  much  more  high-frequency  noise  than 
the  uniform  impulse  in  the  rectangular  membrane  case.  The  noise  above  to*  =  0.25 
occluded  the  response. 


Angular  Frequency  (rad/s) 


(a)  Composite  plot  of  PSDs  for  every  node.  (b)  Comparison  of  model  results  to  exact  so- 

Vertical  lines  depict  exact  modal  frequencies.  lution. 


Figure  15.  Frequency  spectrum  of  the  response  of  the  circular  linear  membrane 
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Response  History  Results. 


The  response  histories  for  four  nodes  of  the  rectangular  membrane  are  shown  in 
Figure  16.  The  thick  line  depicts  the  numerical  prediction.  The  thin  line  depicts 
the  analytical  solution  including  the  first  seven  modes.  The  axis  limits  for  all  of  the 
plots  are  consistent  to  illustrate  relative  magnitudes.  As  expected  from  the  PSD 
in  Figure  14(b),  high-frequency  noise  was  present  in  the  form  of  “jitters.”  Good 
examples  can  be  seen  in  Figure  16(b)  at  four  locations:  the  first  two  peaks,  at  5.5 
seconds,  and  at  12  seconds.  Otherwise  the  model  predictions  closely  matched  the 
exact  solution. 


(a)  Point  (0.1, 0.1)  (b)  Point  (0.3, 0.6) 


(c)  Point  (0.6, 0.5)  (d)  Point  (0.9, 0.8) 

Figure  16.  Response  histories  for  the  rectangular  linear  membrane. 

The  non-uniform  impulse  excited  in  the  circular  membrane  case  was  more  chal- 
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lenging  for  the  numerical  model.  The  responses  of  nodes  at  radius  r  =  0.16,  0.32,  0.54, 
and  0.73  are  depicted  in  Figure  17.  The  thick  line  depicts  the  numerical  prediction. 
The  thin  line  depicts  the  analytical  solution  including  the  first  seven  modes.  The 
axis  limits  for  all  of  the  plots  are  consistent  to  illustrate  relative  magnitudes.  High- 
frequency  noise  was  clearly  evident,  particularly  around  sharp  slope  breaks  of  the 
lower  modes.  Despite  the  high-frequency  oscillations,  the  lower  modes  were  accu¬ 
rately  predicted  throughout  the  run. 
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(a)  r  =  0.16 


Time  Time 

(c)  r  =  0.54  (d)  r  =  0.73 

Figure  17.  Response  histories  for  the  circular  linear  membrane. 


(b)  r  =  0.32 
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Smoothing  the  Response  Histories. 


The  intensive  effort  over  the  decades  to  optimize  integrators  for  dissipating  high 
frequency  content  has  been  absolutely  crucial  for  predictor  algorithms,  where  errors 
due  to  noise  can  accumulate  through  time  steps.  In  this  simultaneous  TCG  method, 
however,  they  did  not  adversely  impact  the  generally  good  performance  of  the  model 
at  the  targeted  frequencies.  The  spurious  modes  had  a  negligible  effect  on  the  energy 
of  the  system  during  the  short  runs  of  this  study. 

So  what  of  the  requirement  that  an  integrator  dissipate  the  unwanted  high- 
frequency  content?  A  previous  TCG  study  demonstrated  that  dissipation  could 
be  handled  outside  of  the  primary  integration  algorithm.  The  desired  effects  were 
achieved  with  a  dual-algorithm  technique  —  an  accurate  high-order  technique  with 
little  to  no  damping  provided  the  bulk  of  the  response  history,  while  a  second  dissi¬ 
pative  lower-order  method  applied  over  several  time  steps  prior  to  a  recorded  event 
stabilized  the  solution  [62],  Depending  on  the  goals  of  the  analysis,  the  present 
method  can  be  employed  as  part  of  a  two-step  process:  accurate  simultaneous  TCG 
solution  followed  by  the  application  of  a  low-pass  filter.  The  low-pass  filter  smooths 
the  response  histories  by  eliminating  the  undesired  high-frequency  content. 

The  options  available  for  a  low-pass  filter  are  endless  and  we  demonstrate  only  one: 
the  Savitzky-Golay  filter  [99,  114,  104],  Used  primarily  in  spectroscopy,  the  Savitzky- 
Golay  filter  minimizes  distortion  of  the  signal  while  reducing  the  noise.  Most  filters 
operate  in  the  frequency  domain;  in  contrast,  the  Savitzky-Golay  filter  operates  in  the 
time  domain  by  performing  local  least-squares  polynomial  fits  throughout  the  time 
series.  Parameters  include  the  order  of  the  polynomial  n  and  the  size  of  the  window 
N  inside  which  the  local  regression  is  performed.  The  window  size,  also  called  the 
filter  length,  is  defined  as  number  of  nodes  in  the  window.  The  window  extends 
symmetrically  ±kAt  from  node  i  where  k  is  an  integer;  hence  N  =  2k  +  1.  Given 


56 


the  minimum  reliable  wavelength  T  from  Eq.  70  and  the  definition  CFL  =  A t/(h/c), 
the  ideal  window  size  in  terms  of  time  increments  is  the  reciprocal  of  the  normalized 
frequency. 


T  _  4tt  (  1  \ 
At  ~  ~3~  \CFL ) 


(73) 


In  converting  the  ideal  window  size  to  an  odd  integer  N  for  use  in  the  algorithm, 
conservative  rounding  procedures  minimized  the  loss  of  desired  content.  The  operator 
floor  returns  the  largest  integer  less  than  the  argument,  and  odd  returns  the  largest 
odd  integer  less  than  the  argument. 


N  =  odd  (floor  (T/At)  —  1) 


(74) 


For  CFL  <  CFLmax  this  technique  ensured  a  minimum  window  size  of  N  —  3.  When 
CFL  <  0.46,  the  large  window  size  resulted  in  excessively  attenuated  peaks.  Table  3 
shows  the  window  sizes  as  a  function  of  CFL.  The  filter  order  was  set  to  one.  It  is 
important  to  note  that  in  general  the  Savitzky-Golay  parameters  can  be  optimized 
to  provide  good  smoothing  of  long-wavelength  signals  or  good  fidelity  of  narrow  peak 
shapes,  but  not  both  simultaneously  [104],  For  this  demonstration  the  former  goal 
was  chosen. 

Table  3.  Savitzky-Golay  window  size  determination. 


CFL  range 

Window  Size  N 

0.46  <  CFL  <  0.6 

7 

0.6  <  CFL  <  0.83 

5 

0.83  <  CFL 

3 

The  effects  of  the  Savitzky-Golay  smoothing  on  the  rectangular  membrane  solution 
are  shown  in  the  frequency  domain  in  Figure  18,  and  in  the  time  domain  in  Figure  19 
for  comparison  to  the  unfiltered  responses  of  Figure  16.  Likewise,  the  smoothed 
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solution  for  the  circular  membrane  is  shown  in  Figure  20  and  in  Figure  21,  which 
may  be  compared  to  the  raw  solution  of  Figure  17.  In  both  cases,  the  first  two  modes 
remained  largely  intact  with  only  slight  attenuation  of  the  peaks.  The  third  and  fourth 
modes  were  slightly  attenuated.  For  the  rectangular  membrane,  the  spurious  jumps 
were  removed  and  the  smoothed  results  very  nearly  matched  the  exact  solutions.  Even 
for  the  more  challenging  circular  membrane,  the  smoothing  significantly  reduced  the 
noise  and  the  lower  modes  tracked  the  exact  solution  closely. 


Figure  18.  Filtered  frequency  spectrum  for  rectangular  membrane. 

The  demonstrated  smoothing  technique  traded  some  accuracy  at  the  peaks  to 
eliminate  most  of  the  noise  while  preserving  the  targeted  frequency  spectra.  The  se¬ 
lection  and  design  of  the  filter  were  guided  by  the  observed  limitations  of  the  model. 
Many  alternative  filtering  techniques  are  available,  but  none  should  be  applied  arbi¬ 
trarily.  We  merely  demonstrated  that  this  particular  solve-then-smooth  process  was 
a  viable  one  for  the  dynamic  cases  in  this  study. 

Summary 

This  study  has  demonstrated  the  effectiveness  of  a  baseline  STCG  formulation 
in  predicting  targeted  modal  frequencies  and  producing  response  histories  of  a  linear 
membrane.  A  methodology  based  on  signal  sampling  concepts  led  to  a  time  increment 
sufficient  for  capturing  frequency  content  in  a  desired  range.  The  time  increment  was 
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Displacement  Displacement 


(a)  Point  (0.1, 0.1)  (b)  Point  (0.3, 0.6) 


(c)  Point  (0.6, 0.5) 


(d)  Point  (0.9, 0.8) 


Figure  19.  Filtered  response  histories  for  the  rectangular 


linear  membrane. 


Figure  20.  Filtered  frequency  spectrum  for  circular  membrane. 
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no  more  restrictive  than  that  of  conventional  methods.  The  spurious  high-frequency 
oscillations  were  consistent  (thus  characterizable)  and  did  not  cause  instability  or  de¬ 
grade  the  lower-frequency  accuracy.  Post-process  filtering  was  not  offered  as  a  panacea 
or  ultimate  replacement  for  an  optimal  algorithm/element  combination.  However,  for 
these  cases,  characterization  of  the  unwanted  frequency  content  led  to  a  smoothing 
methodology  capable  of  removing  much  of  the  noise  without  affecting  the  modal  fre¬ 
quencies.  These  results  were  accomplished  with  fairly  coarse  spatial  meshes  and  the 
simplest  of  elements. 

Although  the  simultaneous  solution  approach  was  analyzed  in  this  work  when 
applied  to  transient  dynamics,  there  is  greater  potential  for  application  to  periodic 
problems.  Example  scenarios  include  periodic  external  forces  or  periodic  boundary 
conditions,  either  of  which  could  conceivably  be  induced  by  a  wing  flapping  mecha¬ 
nism.  In  this  case,  the  simultaneous  solution  approach  turns  an  initial  value  problem 
(IVP)  into  a  boundary  value  problem  (BVP),  enabling  the  employment  of  innumer¬ 
able  finite  element  techniques  in  the  time  domain.  Based  on  the  results  of  this  study, 
the  method  holds  promise  for  such  an  application. 
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IV.  Nonlinear  Analysis  of  Membrane  Dynamics 


For  numerical  simulation  of  nonlinear  membranes,  a  different  approach  was  taken 
than  the  one  that  led  to  the  Simultaneous  Time-Continuous  Galerkin  (STCG)  method 
of  Chapter  III.  There  were  several  rationales  for  this  decision.  First,  the  nonlinear 
model  was  to  include  both  in-plane  and  out-of-plane  displacements,  so  the  size  of 
the  global  matrices  of  the  STCG  would  quickly  exceed  available  computer  hardware 
limitations.  Second,  since  the  system  would  no  longer  be  linear,  timely  convergence 
to  the  global  solution  given  an  arbitrary  initial  guess  became  a  concern.  Lastly,  the 
high-frequency  noise  present  in  the  STCG  solutions  was  not  catastrophic  in  the  linear 
case,  but  could  potentially  become  more  troublesome  for  a  nonlinear  system. 

Instead  of  the  space-time  FEA  approach,  the  more  common  path  was  taken  by 
separately  discretizing  space  and  time  to  employ  a  time-marching  algorithm.  In  this 
way,  novel  discretization  schemes  could  be  devised  separately  for  both  space  and 
time,  independently  studied,  and  finally  combined  to  fulfill  the  ultimate  objectives 
of  this  research.  This  chapter  is  organized  to  follow  this  logic.  First,  the  Hcrmite 
time  interpolation  scheme  is  detailed,  followed  by  the  point  collocation  spatial  dis¬ 
cretization  scheme.  Finally,  the  combination  of  the  two  schemes  to  simulate  dynamic, 
geometrically  nonlinear  membranes  is  examined. 

Hermite  Time  Interpolation 

This  section  describes  a  Hermite  time  interpolation  method  that  differs  from  pre¬ 
vious  related  works  in  the  literature  [109,  81]  by  its  formulation.  Rather  than  using  an 
integral  formulation,  the  proposed  method  determines  the  unique  quintic  polynomial 
interpolation  using  beginning  and  end  point  constraints  on  the  third  derivative  of  the 
function,  called  the  “jerk”  [123].  This  approach  opens  possibilities  for  performance 
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optimization  by  developing  appropriate  constraint  formulas. 

Besides  proposing  an  alternative  formulation,  the  previously-mentioned  works  will 
be  extended  by  calculating  the  precise  stability  boundaries  for  a  linear  unforced  os¬ 
cillator,  exploring  the  fundamental  source  of  the  approach’s  instabilities,  precisely 
calculating  the  dispersion  (relative  period  error),  and  demonstrating  a  procedure  for 
accurately  estimating  local  error.  To  illustrate  implementation  of  the  technique,  a 
fixed-point  iteration  algorithm  will  be  provided  as  a  launching  point  for  general  sys¬ 
tems  of  nonlinear  partial  differential  equations.  Numerical  examples,  including  mul¬ 
tidimensional  and  stiff  systems,  will  demonstrate  the  versatility  of  the  method  and 
the  utility  of  the  local  error  estimation  procedure. 

Time  Discretization. 

Time  is  discretized  by  the  dimensionless  time  r  during  the  time  interval  t°  to 
t1.  The  time  step  size  is  At  =  t1  —  t°.  In  this  portion  of  the  research,  superscripts 
with  a  number  from  0  to  1  indicate  the  discrete  time  at  which  the  value  is  taken  (for 
example,  u°  is  the  value  of  u  at  the  beginning  of  the  time  step,  where  r  =  0). 

T  =  ^e[0,l],i6  [iV1]  (75) 

Consider  the  scalar  equation  L(y)  =  /,  where  the  differential  operator  L  and  forcing 
function  /  may  be  nonlinear,  and  the  solution  y  is  smooth  (y  e  C°°).  The  function 
value  y  and  its  first  two  derivatives  at  the  beginning  and  end  of  the  time  step  are 
collected  into  the  vector  z, 

z  =  {  y(t°)  y(t°)  y(t°)  y{tl)  y( t x)  ijit1)  j  (76) 
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and  re-labeled  as  the  discrete  nodal  variables  displacement  u,  velocity  v,  and  accel¬ 
eration  a. 

z  =  |  u°  v°  a0  u1  v 1  a 1  |  (77) 

The  information  in  z  is  sufficient  to  determine  a  quintic  polynomial  from  t°  to  t 1  of 
the  form  w(t)  =  ^  for  i  =  0 ...  5  with  constant  coefficients  6*  G  3ft.  The  resulting 
six  quintic  Hermite  shape  functions  H(r)  were  derived  in  Ref.  [36]  and  are  shown  in 
Eq.  78.  For  convenience,  the  first,  second,  and  third  derivatives  of  the  Hermite  shape 
functions  are  provided  in  Eq.  79,  80,  and  81,  respectively. 


( 


H(r)  =  { 


v 


1  —  10r3  +  15r4  —  6  r5 
At  (r  —  6r3  +  8r4  —  3r5) 

A*2  (P  -  §r3  +  §r4  - 

10r3  —  15r4  +  6r5 
At  (— 4r3  +  7r4  -  3r5) 
A  t2(ir3-r4  +  ir3) 


\  T 
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(— 30r2  +  60r3 


30r4) /At 


\  t 


H(r)  = 


1  -  18r2  +  32r3  -  15r4 
At  (r  —  | r2  +  6r3  —  |r4) 
(30r2  -  60r3  +  30r4) /At 


> 


—  12r2  +  28 r3  -  15r4 
At(|r2-4r3  +  |r4) 


(78) 


(79) 
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(— 60r  +  180r2  -  120r3) /At2 
(—36 r  +  96r2  —  60r3)  /  At 


T 


H(r)  =  ^ 


1  -  9r  +  18r2  -  10r3 
(60r  -  180r2  +  120r3) /At2 
(—24 r  +  84r2  —  60r3)  /  At 
3 r  -  12r2  +  10r3 


( 


H(r)  =  { 


v 


(-60  +  360t  -  360r2)  /At3 
(-36  +  192r  -  180r2)  /At2 
(-9  +  36r  -  30r2) /At 
(60  -  360r  +  360r2)  /At3 
(-24  +  168r  -  180r2)  /At2 
(3  -  24r  +  30r2)  /At 


\  t 


> 


(80) 


(81) 


The  nodal  approximation  for  the  third  derivative  of  y  is  also  given  a  label  (j  for 
’’jerk”),  so  the  discrete  approximations  of  the  function  y  and  its  first  three  derivatives 
throughout  a  time  increment  t  G  [t0,  ti]  are  then 


y{t)  «  u(t) 

=  H(r)z 

(82) 

y(t )  ~  v(t) 

=  H(t)z 

(83) 

y(t )  ~  a(r) 

=  H(r)z 

(84) 

y(t)  &j{r) 

=  H(r)z 

(85) 

These  interpolation  formulas  define  a  trajectory  with  a  continuous  acceleration  profile. 
The  displacement,  velocity,  and  acceleration  curves  are  intrinsically  constrained  to 
each  other  and  have  been  “pre-integrated”  -  once  the  acceleration  curve  is  uniquely 
determined,  the  velocity  and  displacement  profiles  automatically  follow.  Also,  as  will 
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later  be  demonstrated,  the  polynomials  will  supply  a  wealth  of  information  about  the 
dynamics  occurring  within  the  time  step  that  can  be  exploited  for  error  estimation. 


Constraint  Formulation. 


The  use  of  jerk  constraints  to  define  the  unique  solution  forms  the  foundation  of 
the  present  method.  The  user-defined  formulas  for  jerk  at  times  t°  and  t 1  are  tied  to 
the  Hcrmite  polynomial  interpolation,  thereby  defining  the  beginning  and  end  slopes 
of  the  cubic  acceleration  polynomial.  Thus  a  unique  set  of  interpolation  polynomials 
for  displacement,  velocity,  and  acceleration  are  defined  by  the  Hermite  constraint 
equations 


(86) 


H(0)z  —  j°  =  0 
H(l)z  j1  =0 

It  is  convenient  at  this  point  to  remove  the  acceleration  variable  from  z,  as  in  practice 
a  will  be  calculated  in  tandem  with  j.  The  resulting  matrix  form  of  the  Hermite 
constraint  equations  is 


(87) 


For  some  problems  the  equation  can  be  simplified  by  eliminating  the  a  and  j  terms 
by  substitution  and  solving  the  linear  system  directly.  Otherwise,  moving  the  a 1 
and  j 1  terms  to  the  right-hand  side  provides  an  expression  for  solving  the  system  by 
fixed-point  iteration.  With  no  effect  on  previous  developments,  the  scalar  degrees 
of  freedom  are  now  shown  as  vectors  to  generalize  the  scheme  for  multidimensional 
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problems. 


i 


ir 


1 

0 


7Af2/20  Af3/20 
At/2  At2 /12 


3Af2/20 

At/2 


— At3/30 
-At2/12 


(88) 


The  fixed-point  method  shown  above  can  also  accurately  be  described  as  successive 
substitution  between  two  sets  of  equations  [50],  as  shown  clearly  by  the  following 
algorithm: 

1.  From  u°  and  v°,  calculate  a0  and  j°  (unless  carried  over  from  previous  time 
step) 

2.  Guess  u1  and  v1 

3.  Calculate  a1  and  j1  (discussed  in  the  next  section) 

4.  Until  converged: 

(a)  Calculate  u1  an  v1  using  Eq.  (88) 

(b)  If  converged,  STOP;  else,  update  a1  and  j1 

The  fixed  point  iteration  method  only  requires  the  calculation  of  explicit  vector  formu¬ 
las  -  the  tasks  of  managing  a  Jacobian  and  solving  a  linear  system  at  each  iteration 
is  avoided.  Convergence  will  be  slower  than  the  Newton  method’s  quadratic  rate, 
but  acceleration  methods  are  available  if  required  [135].  Newton  iterations  are  more 
robust  and  may  converge  in  cases  where  fixed  point  iteration  fails.  Generally,  the 
efficacy  of  fixed-point  iteration  versus  Newton  iterations  tends  to  be  problem-specific 
and  comparisons  may  be  found  in  the  literature  [116,  48]. 
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Keeping  with  the  single-step  perspective,  the  initial  guesses  for  u1  and  v1  in  this 
study  were  formed  from  the  available  information  at  time  f°. 


u°  +  v°At  +  ^a°At2  +  |j°At3 

(89) 

v°  +  a0  At  +  If  At2 

(90) 

The  stopping  criterion  was  based  on  a  tolerance  e  for  the  change  in  the  solution  vector 
after  an  update.  Letting  A  indicate  the  change  in  a  value  after  an  iteration,  iterations 
were  stopped  when 

max  (|Au|moa. ,  |  Av|maJ  <  e  (91) 

Since  accuracy  was  the  focus  of  the  study,  tight  tolerances  were  generally  set  without 
consideration  of  the  number  of  function  calls  required. 


Constraint  Definition. 

Given  the  Hermite  interpolation  scheme  just  described,  constraints  will  be  applied 
at  the  end  points  to  obtain  a  unique  solution.  Consider  a  single  time  step  where 
the  initial  displacement  u°  is  known,  and  the  initial  velocity  v°  is  given  or  may  be 
calculated  from  the  first-order  ODE.  Four  unknowns  remain  in  the  vector  z:  a0,  u 1 , 
v1,  and  a1.  Two  equations  may  be  provided  by  constraining  the  accelerations  in 
accordance  with  the  ODE  or  knowledge  of  the  particular  physical  system. 


a0  =  y  (f°,  u°,  v°,  f°) 
a1  =  y  (f1,  u1,  /') 


(92) 


When  posed  as  a  nonlinear  structural  engineering  problem  with  constant  mass  ma¬ 
trix  M,  damping  matrix  C,  stiffness  matrix  K.  and  external  load  F  the  equilibrium 


equation  is 


Ma  +  Cv  +  Ku  +  F  =  0  (93) 

The  familiar  acceleration  constraint  formulas  are  obtained  by  rearranging  the  equi¬ 
librium  equation. 


a0  =  -M'1  (C°v°  +  K°u°  +  F°) 
a1  =  — M^1  (C1v1  +  K1u1  +  F1) 


(94) 


Now  only  u 1  and  v1  are  unknown.  To  obtain  the  final  two  equations,  j°  and  j 1  must 
be  defined.  There  is  great  freedom  in  choosing  jerk  constraints,  as  long  as  desired 
performance  is  achieved.  To  illustrate  this  point,  suppose  requirements  called  for 
an  algorithm  capable  of  unconditional  stability  and  controllable  dissipation.  The 
derivation  could  lead  to  the  jerk  constraints 


?  = 


?  = 


(95) 

(96) 


where  7  and  f3  are  scalar  parameters.  These  jerk  constraints  revert  the  present  for¬ 
mulation  to  the  classic  Newmark  formulation.  (This  correspondence  can  easily  be 
shown  by  side-by-side  comparison  of  the  Newmark  equations  with  Eq.  88.)  For  non- 
dissipative  Newmark  methods,  where  7  =  1/2,  the  constraints  are  identical  and  an 
elegant  symmetry  is  obtained. 


J°  =  J1 


/a1  —  a° 
V  At 


(97) 
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The  Newmark  jerk  constraints  obviously  provide  excellent  results,  but  they  result  in 
artificially  oscillatory  Hcrmite  interpolations  within  the  time  step.  Greater  accuracy 
can  be  achieved  by  defining  jerk  according  to  its  physical  and  mathematical  definition: 
as  the  time  derivative  of  the  acceleration  constraints  in  Eq.  (92). 


j°  =  V  (t° ,  u°,  v°,  a0,  f°,  /°j 

j1  =  V  (V^'tdV,/1/1) 


(98) 


The  jerk  constraints  are  stated  here  in  terms  of  acceleration,  which  may  be  eliminated 
by  substitution  of  Eqs.  (94)  if  desired. 


j°  =  —AT1 

C°a°  +  (c°  +  K°)  v°  +  K°u°  +  F° 

j1  =  -M-1 

CV  +  (C1  +  K1)  v1  +  K1!!1  +  F1 

(99) 


For  the  remainder  of  this  dissertation,  although  the  Hermite  interpolation  method 
has  been  shown  to  be  a  family  of  methods  depending  on  the  chosen  jerk  constraints, 
the  method  will  be  confined  to  the  jerk  constraints  of  Eqs.  (98). 


Stability  and  Energy  Conservation. 

The  linear,  undamped,  free  system  y  +  c o2y  =  0  with  constant  angular  frequency 
oj  provides  the  prototypical  case  for  investigating  the  method’s  stability  and  energy 
conservation.  The  constraint  equations  trivially  fall  from  the  ODE  as  a  =  — c u2u 
and  j  =  —u2v.  The  constraint  equations  are  substituted  into  Eq.  88  to  produce  the 
iteration  matrix  A,  which  for  a  time  step  maps  the  initial  state  {u°,  v°}T  to  the  final 
state  {u1,v1}T.  The  time  step  size  is  non-dimensionalized  as  v  =  ujAt  =  27rAf/T, 
where  T  is  the  period  of  the  oscillator.  While  v  simplifies  the  matrices,  most  of 
the  following  discussion  will  reference  A t/T,  due  to  its  widespread  use  in  the  finite 
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element  literature.  The  iteration  matrix  simplifies  to 


A 


1 

z/4  +  16z/2  +  240 


3z/4  -  104z/2  +  240  (J)  (|z/5  -  24z/3  +  240z/) 

-uj  (— 24z/3  +  240z/)  3z/4  -  104z/2  +  240 


(100) 


which  is  the  method’s  approximation  of  the  exact  iteration  matrix 


A 


exact 


cos  H  (£) sin  H 

— ca  sin  (i/)  cos  (u) 


(101) 


The  exact  iteration  matrix  is  typically  found  in  Hamiltonian  formulations  using  the 
conjugate  momentum  p,  which  is  equivalent  to  this  formulation  by  assuming  without 
loss  of  generality  that  the  mass  m  —  1,  so  v  —  p.  [116,  48].  For  this  unforced  linear 
oscillator,  the  present  method  produces  an  iteration  matrix  with  the  same  diagonal 
elements  as  the  weighted  residual  method  of  Razavi  et  al.  [109].  Therefore,  the  sta¬ 
bility  characteristics  will  be  the  same  and  their  findings  can  be  further  developed.  An 
in-depth  examination  of  the  iteration  matrix  follows,  to  provide  greater  insight  into 
the  method. 

The  components  of  the  exact  iteration  matrix  and  its  approximation  are  compared 
in  Figure  22.  The  horizontal  axes  represent  a  chosen  step  size  A t/T,  and  the  vertical 
axis  is  the  value  of  the  corresponding  matrix  element.  The  dashed  lines  are  used  for 
the  exact  matrix,  Aexact,  and  the  solid  lines  are  used  for  the  approximate  matrix, 
A.  The  plus  sign  markers  show  the  values  of  the  iteration  matrix  obtained  from 
Ref.  [109].  The  absolute  value  of  the  errors  is  shown  in  the  bottom  right  subplot 
up  to  the  Nyquist  frequency  (note  the  change  to  a  logarithmic  vertical  axis).  The 
approximations  are  quite  accurate  for  reasonable  time  steps,  with  error  magnitudes 
of  the  matrix  elements  less  than  1  x  10~4  for  A  t/T  =  0.1.  However,  the  downside 
of  using  this  polynomial  approximation  is  also  obvious,  as  with  increasing  time  step 
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size,  the  approximations  exceed  the  bounding  range  of  the  exact  values  and  A12  grows 
linearly  without  bound.  The  regions  of  instability  defined  by  this  behavior  will  be 
characterized  next. 


Figure  22.  Graphical  depiction  of  the  Hermite  time  interpolation  method  iteration 
matrix  for  the  simple  harmonic  oscillator 


A  symplectic  method  accurately  reproduces  the  geometric  structure  of  the  ODE 
and  its  solutions,  and  conserves  invariants  over  long-duration  simulations  [48].  For 
Hamiltonian  systems  in  mechanics  the  invariant  is  typically  energy.  It  is  easily  proven 
that  the  iteration  matrix  A  is  symplectic  because  for  all  u>  and  v  it  satisfies  the 
relationship 

AtJA  =  J  (102) 

where  J,  the  skew-symmetric  operator  used  when  representing  a  Hamiltonian  system 


72 


of  equations  in  matrix  form,  is  defined  as 


J 


0  I 
-/  0 


(103) 


As  with  any  symplectic  matrix,  it  is  also  true  that  det(A)  =  1  [89].  Since  A  is 
an  area-preserving  linear  mapping,  the  stability  can  be  determined  by  finding  where 
[hr (A)  |  <  2  [2];  or,  since  An  =  A2 2,  where  |An|  <  1.  Hence  the  stability  can  be 
inspected  in  the  upper-left  plot  of  Figure  22.  It  is  immediately  observed  that  there 
are  two  regions  of  instability  -  a  small  region  just  above  A t/T  =  0.5,  and  a  clear  limit 
above  A  t/T  «  1.2,  corresponding  to  those  of  Ref.  [109].  The  regions  of  instability 
are  shaded  gray  in  the  figure.  We  now  determine  the  exact  locations  of  the  stability 
limits.  The  boundaries  are  determined  from  the  polynomial  inequality 


3z/4  -  104z/2  +  240 


(104) 


The  terms  with  v  can  be  lumped  into  a  single  term, 

2z/4  -  120z/2 
+  ~a  +  I61/2  +  240 


(105) 


with  some  final  manipulations  showing  that  the  stability  requires  satisfaction  of  two 
inequalities. 

2z/4  -  120u2 
z/4  +  16z/2  +  240 


z/4  —  60z/2 
z/4  +  16z/2  +  240 


<  0 


(107) 


The  first  statement  is  violated  if  10  <  v2  <  12  and  the  second  is  violated  if  v2  >  60, 
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providing  the  exact  instability  regions  in  terms  of  sampling  frequency  as 


V5/(2t t2)  <  A t/T  <  y/3 /tt2  (108) 

V15/tt2  <  At/T  (109) 


The  instability  can  be  further  characterized  by  an  eigenvalue  analysis.  The  expression 
for  the  eigenvalues  is 

_  3z/4  -  104z/2  +  240  ±  2z/a/2  (z/2  -  60)  (z/2  -  12)  (z/2  -  10) 

*  “  z/4  +  16z/2  +  240 


and  the  stability  boundaries  can  be  seen  by  inspection  of  the  expression  under  the 
radical.  When  the  radicand  is  positive  the  eigenvalues  lie  on  the  real  axis  of  the 
complex  plane  and  instability  results  [28].  The  path  of  the  eigenvalues  as  the  time 
step  is  increased  is  depicted  in  Figure  23(a).  Starting  at  the  dots  on  the  right- 
hand  side  and  increasing  the  step  size,  the  eigenvalues  follow  the  unit  circle  in  the 
negative-real  direction  until  they  collide  and  separate  at  (—1,0),  corresponding  to 
A  t/T  =  \/5/(27r2).  They  quickly  return  to  the  circle  and  follow  it  to  (1,0),  at 
which  point  they  again  collide  and  repel  each  other.  They  remain  on  the  real  axis  as 
A  t/T  — *  oo.  For  a  detailed  discussion  of  this  behavior,  see  Ref.  [2], 

The  moduli  of  the  two  eigenvalues  as  a  function  of  A  t/T  are  plotted  in  Fig¬ 
ure  23(b).  The  spectral  radius  p( A),  defined  as  the  maximum  of  the  moduli,  is  equal 
to  one  wherever  the  method  is  stable  -  no  numerical  damping  is  present  and  the 
energy  is  conserved.  The  traces  also  show  that  after  bifurcation,  the  moduli  are  the 
reciprocals  of  each  other,  another  indicator  of  a  symplectic  transformation  [2],  Ex¬ 
ploring  the  nature  of  the  instability  provides  a  launching  point  for  devising  mitigation 
techniques.  A  small  amount  of  physical  damping  can  eliminate  the  small  unstable 
region  just  above  the  Nyquist  frequency  [109].  However,  physical  damping  alone  af- 
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(a)  Root  locus  trace  of  the  two  eigenvalues.  (b)  Spectral  radius  of  the  iteration  matrix  A. 

Figure  23.  Eigenvalue  behavior  for  the  simple  harmonic  oscillator. 

fects  only  middle  frequencies  and  cannot  remove  high  frequency  instabilities  [58]  - 
algorithmic  damping  would  be  required,  a  topic  left  for  future  research. 


Dispersion. 

To  quantify  the  method’s  dispersion,  let  T  be  the  period  of  the  model  solution  and 
T  be  the  exact  period.  Then  the  period  error  is  P  =  T /T  [28]  and  the  relative  period 
error  is  (T  —  T)/T  [58].  The  period  error  P  can  be  calculated  from  the  eigenvalues 
of  Eq.  110  as 


„  T 
P  =  -  =  v 
T 


'  /  W2^-60)^-12)(y»-10)- 

\  3z/4  -  104z/2  +  240 


-i 


(111) 


As  shown  in  Figure  24(a),  the  Hermite  method  accurately  reproduces  the  frequency 
of  the  oscillator  up  to  the  Nyquist  frequency  with  a  maximum  relative  period  error  of 
0.030.  The  error  is  less  than  0.0016  below  A t/T  =  0.2.  The  central  difference  shrinks 
the  period  until  it  hits  its  stability  limit  at  A  t/T  =  1/n  [28].  The  average  accelera¬ 
tion  method  significantly  increases  the  period.  To  illustrate  the  extent,  Figure  24(b) 
portrays  the  period  errors  for  a  simple  harmonic  oscillator  using  the  Hermite  and 
average  acceleration  methods  at  A  t/T  =  0.4.  Because  the  numerical  samples  are  so 
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sparse,  least-square  fits  of  sinusoidal  functions  to  each  method  are  included  in  the 
figure.  The  fidelity  of  frequency  reproduction  bodes  well  for  the  application  of  this 
method  in  a  wave  propagation  problem. 


At/T 

(a)  Relative  period  errors. 


(b)  Example  of  period  error  at  At/T  =  0.4. 


Figure  24.  Relative  period  error  for  the  simple  harmonic  oscillator  comparing  various 
integration  methods. 


Local  Error  Estimation. 

Error  estimates  can  alert  the  user  to  potentially  inaccurate  solutions,  trigger 
damping  mechanisms,  or  direct  step  size  changes.  A  common  technique  is  to  compare 
solutions  from  two  numerical  methods  of  different  order  [3].  The  Runge-Kutta  4/5 
has  proven  useful  because  it  does  so  with  only  one  more  function  call  [117].  Since 
the  present  method  enforces  constraints  at  the  end  points  of  a  time  step  to  dictate 
the  solution,  the  degree  to  which  those  constraints  are  violated  in  the  interior  offers 
useful  data  for  error  estimation.  An  example  will  be  developed  here  which  uses  one 
additional  function  call  at  mid-step  to  estimate  the  local  error. 

Suppose  the  solution  for  a  time  step  has  been  obtained  so  the  vector  z  is  known. 
Let  the  interpolated  value  of  acceleration  at  the  midpoint  r  =  0.5  be  a0'5  =  H(0.5)z. 
Also  at  the  midpoint,  the  acceleration  constraint  equation  is  applied  to  compute  the 
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constraint  value  of  acceleration  a0'5  =  y  (t0  5,  u0'5,  v0'5,  /°'5).  Thus  the  acceleration 
error  8a  is  the  difference  between  the  constrained  acceleration  and  the  interpolated 
acceleration  at  the  midpoint  of  the  time  step. 

8a  =  a0'5  -  aa5  =  H(0.5)z  -  y  (ta5,  u a5,  ua5,  /°'5)  (112) 

For  a  consistent  formulation,  as  the  time  step  approaches  zero,  the  interpolation 
polynomials  approach  the  exact  solution  throughout  the  time  step  and  8a  approaches 
zero.  The  acceleration  error  will  now  be  used  to  build  an  alternative  acceleration 
profile  which  will  be  integrated  twice  to  produce  a  useful  error  parameter  in  terms 
of  displacement.  First,  three-point  Lobatto  integration  puts  the  acceleration  error 
in  terms  of  velocity.  The  alternative  acceleration  polynomial  is  assumed  to  pass 
through  the  the  points  (t°,a°),  (t°"5,a05),  and  (t1,^1)  where  a1  =  a1  +  8a.  The 
Lobatto  integration  formula  then  provides  a  new  estimate  for  the  final  velocity,  v\. 

vl  =  v°  +  y  ^a° +  ^°  5  +  5fil)  (113) 

Finally,  the  local  error  estimate  for  displacement,  8u.  is  calculated  by  a  simple  ap¬ 
proximation  of  the  integral  of  the  velocity  error  through  the  time  step. 

(114) 

To  correlate  the  local  error  estimate  with  the  exact  local  error,  numerical  experiments 
were  conducted  using  a  simple  harmonic  oscillator.  For  the  experiment,  1,500  data 
points  were  created  by  sampling  each  of  the  time  step  sizes  A t/T  =  0.05,  0.1,  0.2,  0.3, 
and  0.4  with  300  samples  for  the  phase  (j)  =  [0,  2n\  where  y  =  cos(t  —  4>).  A  scatter 
plot  of  the  relationship  is  shown  in  Figure  25.  Nearly  all  of  the  1,500  data  points  he 
near  the  line  of  slope  one.  As  the  error  estimate  8u  increases  past  0.001,  the  error  is 


8u  = 


At  , 
■5r( V 


>i) 
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Figure  25.  Relationship  between  the  calculated  local  error  estimate  and  the  exact  local 
error  for  the  simple  harmonic  oscillator. 

increasingly  overestimated.  The  sets  of  points  for  each  A t/T  have  a  forked  appearance 
from  the  process  of  converting  elliptical  plots  of  the  estimated  versus  exact  error  to 
the  log-log  plot  shown  here.  The  error  estimation  process  just  described  was  used  for 
the  numerical  examples  of  this  paper  because  of  the  strong  correlation  between  the 
estimated  and  exact  local  error  for  reasonable  time  step  sizes  and  the  conservative 
overestimation  of  error  at  extreme  step  sizes. 

Numerical  Examples. 

Simple  Harmonic  Oscillator. 

The  simple  harmonic  oscillator  discussed  earlier  was  solved  from  two  starting  con¬ 
ditions:  (2/(0) ,  2/(0))  =  (1,0)  and  (0,1).  For  comparison,  solutions  were  also  obtained 
using  the  two  members  of  the  Newmark  family:  the  central  difference  method,  and 
the  trapezoid  method  with  7  =  1/2  and  /3  —  1/4  [102],  Convergence  for  all  three 
methods  is  shown  in  Figure  26,  and  the  Hermite  interpolation  accuracy  is  clearly 
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superior  to  the  second-order  Newmark  methods.  The  observed  rate  of  convergence 
matched  the  rates  observed  using  other  Hermite  polynomial  methods  [109,  81]. 


(a)  y( 0)  =  1,  2/(0)  =  0  (b)  2/(0)  =  0,  2/(0)  =  1 

Figure  26.  Convergence  for  the  simple  harmonic  oscillator  from  two  different  initial 
conditions. 


As  expected  after  finding  the  iteration  matrix  to  be  symplectic  by  way  of  Eq.  (102), 
the  energy  of  the  oscillator  was  conserved  over  long  durations.  Figure  27  shows  the 
energy  error  at  the  end  of  1000  periods,  with  the  left  plot  depicting  the  last  ten 
periods  with  A t/T  =  0.1,  and  the  right  plot  depicting  the  last  100  periods  with 
A t/T  =  0.3.  The  error  oscillates  but  remains  bounded  in  a  tight  range,  typical 
behavior  for  symplectic  methods  [68]. 

Nonlinear  Second-Order  ODE. 

In  Ref.  [3],  the  rates  of  convergence  of  several  techniques  were  compared  by  solving 
a  second-order  nonlinear  ODE.  The  present  method  will  now  be  applied  to  the  same 
problem  with  two  objectives:  (1)  confirm  the  order  of  convergence  for  the  method, 
and  (2)  demonstrate  its  application  in  a  fixed-point  iterative  algorithm.  The  ODE  is 


V  =  —e 


y+ 1 


(115) 
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(a)  A t/T  =  0.1  (b)  A t/T  =  0.3 

Figure  27.  Energy  errors  at  the  conclusion  of  a  1000-period  simulation  for  two  different 
time  step  sizes. 


with  the  initial  conditions  y( 0)  =  0  and  y( 0)  =  #tanh(#/4).  By  setting  the  parameter 
6  =  3.03623184819656,  the  exact  solution  for  the  convergence  study  is  y(  1)  =  0.  The 
ODE  directly  provides  the  acceleration  constraints  of  Eq.  (92).  The  time  derivative 
of  the  ODE  provides  the  jerk  constraints  of  Eq.  (98). 


a 

3 


=  — e“+1 


=  —  ve 


u-\- 1 


(116) 

(117) 


These  constraint  formulas  contribute  the  acceleration  and  jerk  values  to  the  right- 
hand  side  of  Eq.  88.  The  resulting  matrix  equation,  shown  in  expanded  form  in 
Eq.  118,  was  used  in  the  fixed-point  iteration  algorithm  to  obtain  the  solutions.  The 
subscript  k  indicates  current  values,  while  k  +  1  indicates  the  new  values  upon  com¬ 
pletion  of  the  iteration.  The  errors  and  rates  of  convergence  in  Table  4  confirm  the 
fourth-order  convergence  seen  in  the  simple  harmonic  oscillator.  In  comparison  to  the 
methods  presented  in  Ref.  [3],  the  error  magnitudes  were  approximately  the  same  as 
those  of  the  fourth-order  Runge-Kutta  method,  but  not  quite  as  good  as  the  four-step, 
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fifth-order  Adams-Moulton  method. 


60 

At3 

60 

At3 


(118) 


Table  4.  Convergence  results  for  the  exponential  function. 


log  (At) 

log(error) 

P 

0 

-1.01 

-1 

-5.28 

4.27 

-2 

-9.27 

3.99 

-3 

-13.26 

3.99 

System  of  Nonlinear  Equations. 

The  next  example  was  chosen  to  demonstrate  the  present  method’s  applicability 
to  systems  of  equations  and  compare  the  error  estimation  accuracy  to  other  tech¬ 
niques.  This  particular  problem  was  taken  from  Reference  [77],  where  it  is  used  to 
detail  accuracy  and  error  estimation  for  several  integration  techniques.  For  direct 
comparison  to  the  author’s  presentation,  the  local  error  was  calculated  as  follows. 
For  each  time  step  (not  just  the  first  one),  the  initial  values  u®  and  v®  were  set  to  the 
exact  solution.  Thus,  at  the  conclusion  of  the  time  step,  the  exact  local  error  was 
precisely  the  difference  between  the  model’s  solution  (uj  and  vj)  and  the  exact  solu¬ 
tion  (y^t1)  and  ijiit1)).  The  exact  local  error  can  then  be  compared  to  the  estimate 
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provided  by  the  model.  The  stated  error  is  the  L2  norm  of  the  displacement  error 
vector.  The  system  consists  of  two  nonlinear  equations,  with  t  G  [0,6.4]: 


in  =  -yi  (: vl  +  yl)  3/2 

(119) 

in  =  -7/2  (7/1  +  t/2)  3/2 

(120) 

The  initial  conditions  7/1  (0)  =  1,  7/1  (0)  =  0,  7/2(0)  =  0,  and  t/2(0)  =  1  provide  the  exact 
solution  7/1  =  cos (t)  and  y2  =  sin(f).  As  with  the  previous  example,  the  ODE  directly 
provides  the  acceleration  constraints  and  the  time  derivative  of  the  ODE  provides  the 
jerk  constraints  to  be  applied  at  the  beginning  and  end  of  the  time  steps  (Eq.  (92) 
and  Eq.  (98),  respectively).  The  superscript  time  indices  are  omitted  for  clarity. 


Oi  =  — 77i  (u\  +  7/2)  3//2  (121) 

a2  =  -772  (ul  +  ul) 3/2  (122) 

jl  =  77 1  (ul  +  uj)  3/2  +  5771  (77177!  +  U2V2)  (l u\  +  77^)  (123) 

32  =  -V2  (t7i  +  772)  3/^  -  772  (t7i77i  +  U2V2)  (u\  +  77^)  ^  (124) 


The  system  was  solved  using  the  fixed-point  iteration  scheme  (Eq.  88)  with  an  itera¬ 
tion  stopping  tolerance  of  1  x  10”10.  The  results  are  shown  in  Table  5.  The  present 
method  was  more  accurate  than  the  reference’s  examples  by  an  order  of  magnitude 
or  more.  The  local  error  estimates  underestimated  the  exact  local  error  by  an  average 
magnitude  of  7%,  5%,  and  6%  for  At  =  0.8,  0.4,  and  0.2  respectively.  By  compar¬ 
ison,  for  At  =  0.2,  the  Fchlberg  (4,5)  Runge-Kutta  method  average  local  error  was 
3.7  x  10~6  and  the  error  estimate  was  off  by  an  average  of  18%  [77]. 

To  check  convergence  rates  of  the  method  once  again,  the  simulation  was  next 
run  without  performing  the  exact  solution  reset  used  above  for  local  error  estimation 
analysis.  The  errors  were  thus  the  difference  between  the  model  results  and  the 
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Table  5.  Local  error  E  and  local  error  estimates  Eest  for  the  system  of  ODEs. 


At  = 

0.8 

At  = 

0.4 

At  = 

0.2 

t 

E 

Eest 

E 

Eest 

E 

Eest 

X 

O 

xlO4 

xlO6 

xlO6 

xlO8 

xlO8 

0.8 

1.83 

1.65 

2.79 

2.61 

4.41 

4.14 

1.6 

1.44 

1.52 

2.64 

2.57 

4.36 

4.13 

2.4 

1.48 

1.53 

2.72 

2.59 

4.40 

4.14 

3.2 

1.83 

1.63 

2.85 

2.63 

4.45 

4.15 

4.0 

1.81 

1.64 

2.77 

2.61 

4.41 

4.14 

4.8 

1.42 

1.51 

2.64 

2.57 

4.36 

4.13 

5.6 

1.50 

1.53 

2.73 

2.59 

4.40 

4.14 

6.4 

1.85 

1.64 

2.86 

2.63 

4.45 

4.15 

exact  solution  at  the  final  time,  t  =  6.4.  Convergence  rate  calculations  are  shown 
in  Table  6.  The  observed  rates  of  convergence  approached  four  as  the  time  steps 
decreased,  confirming  the  findings  from  the  previous  examples. 

Table  6.  Convergence  results  for  the  multidimensional  example.  The  order  of  conver¬ 
gence  is  labeled  p. 


At 

log(At) 

log(y1  error) 

P 

log  (1/2  error) 

P 

log(L2  error) 

P 

0.8 

-0.097 

-2.40 

- 

-2.76 

- 

-2.36 

- 

0.4 

-0.398 

-3.54 

3.76 

-3.86 

3.67 

-3.49 

3.74 

0.2 

-0.699 

-4.72 

3.95 

-5.04 

3.92 

-4.68 

3.94 

Stiff  System  with  Variable  Step  Sizes. 

For  the  final  example,  a  stiff  system  was  selected  to  challenge  the  error  estimation 
procedure  and  demonstrate  its  effectiveness  for  step  size  control.  Briefly  stated,  a 
stiff  system  typically  requires  exceedingly  small  time  steps  to  resolve  with  a  stability- 
limited  explicit  method,  even  though  the  solution  may  appear  smooth.  Often,  both 
slow  and  fast  responses  (transients)  are  present.  Fixed-point  iteration  typically  strug¬ 
gles  to  converge  [77].  Further  explanation  may  be  obtained  from  Refs.  [3]  and  [49], 
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which  both  discuss  the  example  selected  here.  It  features  a  rapid  initial  transient 
followed  by  a  slower  response  that  loosely  follows  y  =  cos(t).  The  first-order  ODE  is 

y  =  — 100  (y  —  cos(f)) ,  f  e  [0,  7t/2]  (125) 

from  which  the  acceleration  and  jerk  constraints  can  be  formulated  in  terms  of  v  and 
t. 


a  =  — lOOu  —  100sin(f)  (126) 

j  =  lOOOOu  +  10000  sin(f)  -  100  cos (t)  (127) 


The  fixed-point  iteration  algorithm  used  so  far  in  this  paper  required  the  domain 
to  be  divided  into  at  least  100  time  steps  to  converge,  and  the  accuracy  was  still 
poor.  Fortunately,  the  system  of  equations  is  linear  with  respect  to  velocity  and  the 
acceleration  and  jerk  terms  can  be  eliminated  by  substitution.  The  time-dependent, 
nonlinear  forcing  terms  are  lumped  into  the  vector  b,  resulting  in  the  system 
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+  (b(t°,t1)} 


where  the  forcing  vector  is 


(128) 


(129) 


n  ,  [  (-w7  +  10000)sin(f0)-100cos(f°)  +  ^sin(f1)  I 

b{t°,t1)  =  {  y  At  ’  At  }  (130) 

I  —  ^sin(f°)  +  +  10000)  sin(f1)  —  100cos(f1)  j 
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For  each  new  time  step,  t°,  t1,  and  At  are  known;  the  coefficient  matrices  can  be 
calculated  once  per  simulation  (or  per  time  step  if  At  is  varied),  and  the  vector 
b(t°,t1)  can  be  calculated  once  per  time  step.  The  resulting  linear  system  can  be 
solved  directly. 

First,  the  problem  was  solved  using  fifty  uniform  time  steps  of  At  =  7r/ 100. 
The  results  are  shown  in  Figure  28(a).  The  dots  are  nodal  values  of  the  present 
method.  The  solid  red  line  is  the  benchmark,  a  converged  solution  produced  by 
the  Scientific  Python  (SciPy)  function  “scipy . integrate . odeint”  [65].  The  SciPy 
function  uses  variable-order  Adams  and  backward  difference  formula  (BDF)  routines 
from  the  FORTRAN  library  “odepack"  and  can  accommodate  stiff  systems.  The 
cumulative  local  error  estimates  are  displayed  as  bands  above  and  below  the  present 
model’s  solution. 

The  model  clearly  overshoots  the  initial  level-off  before  following  the  ODE  con¬ 
tours  parallel  to  the  benchmark  solution.  By  itself  the  model’s  solution  is  misleading 
-  it  is  qualitatively  correct  but  quantitatively  inaccurate.  Using  the  benchmark  as 
the  truth,  the  global  error  at  the  last  point  was  0.135. 

The  error  estimation  process  successfully  detected  the  struggle  of  the  model  in 
following  the  rapid  transient  and  widened  the  error  band  appropriately.  As  the  solu¬ 
tion  settled,  the  local  error  estimates  decreased  significantly  (in  other  words,  the  error 
band  width  did  not  shrink  or  expand  as  time  progressed).  The  benchmark  solution 
happens  to  lie  within  the  cumulative  local  error  bands,  though  this  may  not  always 
be  the  case,  and  the  cumulative  local  error  bands  should  not  be  construed  as  global 
error  bounds.  In  the  absence  of  the  reference  solution,  the  wide  error  band  would 
suggest  to  the  user  that  the  solution  is  not  trustworthy. 

Besides  suggesting  the  solution  is  inaccurate,  the  error  estimates  can  also  be  used 
to  vary  the  step  sizes.  To  demonstrate  this  utility,  the  problem  was  re-run  with  step 
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size  control  performed  by  a  simple  rule  from  Ref.  [117].  Let  e  be  the  user-defined  upper 
bound  for  the  local  truncation  error.  The  local  error  estimate  Su  is  obtained  from 
Eq.  (114)  and  measured  using  the  maximum  norm.  After  a  time  step  is  completed,  a 
new  step  size  is  calculated  using  the  formula 


At 


new 


IMI 


1/5 


At 


previous 


(131) 


If  the  tolerance  for  the  current  step  is  violated,  the  new  step  size  is  applied  to  the 
current  step  and  it  is  re-calculated.  If  the  tolerance  is  not  violated,  the  current  step 
solution  is  accepted  and  the  new  step  size  is  applied  to  the  next  step. 

The  solution  with  variable  step  sizes  is  shown  in  Figure  28(b).  An  error  tolerance 
of  0.01  was  used  for  accepting  a  time  step’s  solution;  fairly  loose,  but  effective  for 
displaying  the  new  error  estimate  lines  on  the  plot.  The  model  repeated  a  step  only 
once  (the  first  time  step)  and  appropriately  grouped  the  points  at  the  initial  peak. 
Only  twelve  steps  were  performed  and  the  error  was  reduced  from  0.135  to  0.0175. 
Again,  the  benchmark  happens  to  he  within  the  representative  error  bands.  More 
significantly,  the  reported  errors  are  significantly  smaller  than  in  the  first  example, 
signalling  a  more  accurate  solution. 

A  final  case  was  run  with  a  step  size  error  tolerance  of  10”8.  The  last  point’s  global 
error  was  2.17  x  10-6  with  about  the  same  number  of  steps  (52)  as  the  constant  step 
size  run  (50).  On  four  occasions  during  the  simulation,  the  algorithm  estimated  the 
local  error  to  be  out  of  tolerance,  and  therefore  the  time  step  was  repeated  with 
a  smaller  step  size.  The  local  error  estimate  lines  were  indistinguishable  from  the 
benchmark  solution  so  no  plot  is  shown. 


(a)  Constant  time  steps  (b)  Variable  time  steps 

Figure  28.  Stiff  system:  Solution  and  error  bands  compared  to  the  benchmark. 

Summary. 

This  section  has  demonstrated  the  effectiveness  of  solving  initial  value  problems 
with  quintic  Hermite  polynomial  interpolations  defined  by  end-point  jerk  constraints. 
The  Newmark  methods  were  shown  to  be  a  subset  of  the  current  formulation.  Com¬ 
plete  analysis  was  performed  with  the  jerk  constraints  defined  according  to  the  gov¬ 
erning  equation.  The  method  is  superior  to  the  second-order  Newmark  methods  in 
terms  of  absolute  accuracy,  rate  of  convergence,  and  frequency  reproduction.  The 
conditional  stability  was  fully  characterized  using  a  linear  oscillator,  and  it  was  found 
that  the  regions  of  instability  existed  only  above  the  Nyquist  frequency  -  well  above 
the  time  step  sizes  demanded  by  accuracy  requirements.  Nonlinear  systems  are  readily 
accommodated.  The  problem  formulation  is  systematic  and  physically  intuitive,  and 
can  be  efficiently  executed  in  a  fixed-point  iteration  algorithm  for  non-stiff  systems. 

The  accuracy  and  application  of  the  method  has  been  emphasized,  not  the  compu¬ 
tational  efficiency.  During  the  course  of  the  study  and  in  initial  testing  of  the  nonlinear 
membrane  model,  the  Hermite  interpolation  method  solved  by  the  fixed-point  itera¬ 
tion  algorithm  has  been  suitably  efficient.  Considering  the  accuracy  attained  versus 
the  computational  expense  (both  function  calls  and  run  time),  it  is  believed  that  this 
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method  can  be  competitive  with  other  time  integration  techniques. 


Point  Collocation  Spatial  Discretization 

In  this  section,  a  novel  membrane  model  based  on  the  group  finite  element  for¬ 
mulation  and  point  collocation  method  will  be  presented  and  evaluated  for  static 
membranes.  After  summarizing  the  governing  equations  and  the  material  model, 
polygon  interpolation  formulas  will  be  derived  to  calculate  gradients  in  a  staggered 
grid  approach.  The  steps  for  calculating  nodal  force  imbalances  will  be  described  in 
detail.  After  the  model  is  explained,  verification  will  demonstrate  consistency  and 
an  observed  rate  of  convergence  of  two.  Finally,  predictions  will  be  validated  against 
experimental  results  in  the  literature  to  show  the  model  to  be  suitable  through  its 
range  of  intended  use  (he.,  short  of  the  onset  of  hyperelastic  material  response). 


Governing  Equations. 

For  this  study,  a  membrane  is  defined  as  a  thin  plate  without  bending  stiffness  [63] . 
The  governing  partial  differential  equations  (PDEs)  for  the  nonlinear  membrane  are 
those  of  a  plate  undergoing  finite  deformations  as  derived  in  [111],  but  with  the  mo¬ 
ment  and  curvature  terms  removed.  The  internal  stress  components  are  Nij  (second 
Piola-Kirchoff,  per  unit  length)  with  the  subscripts  denoting  orientation  with  respect 
to  the  Cartesian  axis  directions  X\  and  x2.  The  external  pressure  vector  components 
are  represented  by  /.  The  terms  are  expressed  in  force  per  unit  area. 

(132) 

(133) 
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The  solution  of  the  membrane  surface  displacements  u  begins  by  recovery  of  the 
surface  gradients.  The  gradients  then  lead  to  the  components  of  the  Green  strain 
tensor  E \j  [111,  15]. 
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The  internal  stresses  are  calculated  using  the  conventional  plane  stress  constitutive 
relationship.  Prestress  is  accounted  for  by  the  vector  No-  Note  that  the  membrane 
thickness  h,  where  h  =  h( u),  is  a  function  of  the  displacement  field  to  take  thinning 
into  account  as  the  membrane  stretches. 
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Discretization  and  Polygon  Interpolation  Formulas. 

By  using  the  group  formulation,  all  three  governing  PDEs  (Eq.(132),  Eq.(133), 
and  Eq.(134))  were  cast  into  the  same  first-order  PDE  form 
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where  the  vectors  Q  define  the  degrees  of  freedom  as 


The  domain  is  discretized  by  forming  a  staggered  mesh,  meaning  different  variables 
are  evaluated  at  different  points  in  the  domain.  For  example,  the  displacements  u  and 
the  stresses  N  will  not  be  computed  at  the  same  nodes.  As  commonly  found  in  finite 
difference  discretizations,  a  staggered  mesh  enables  more  compact  stencils.  In  certain 
CFD  applications,  high-frequency  oscillations  are  reduced  because  the  pressure  and 
velocity  fields  are  fully  coupled  [55]. 

The  staggered  mesh  consists  of  a  primary  mesh  of  three-node  linear  triangles  and  a 
dual  mesh  of  polygons.  The  nodes  of  the  primary  mesh  define  the  model’s  collocation 
points  and  carry  the  vectors  u  (displacement)  and  f  (external  force).  The  role  of 
the  primary  mesh  is  to  recover  the  first  partial  derivatives  of  the  membrane  surface. 
The  calculated  partials  of  each  triangle  are  placed  at  the  centroid,  as  is  common  in 
post-processing  gradient  recovery  procedures  [28,  146,  98]. 

The  dual  mesh  is  formed  by  connecting  the  centroids  of  the  triangles  to  form 
vertex-centered  polygons,  also  called  tributary  areas  [146].  Note  that  the  polygons  do 
not  overlap  and  are  therefore  not  the  same  as  an  element  patch.  The  vertices  of  the 
polygons  carry  the  Q  vectors.  The  solution  of  the  discretized  governing  equations, 
Eq.  (137),  requires  an  approximation  of  the  partial  derivatives  of  Q  at  the  center 
node.  For  this  study,  the  polygon  patch  interpolation  presented  in  [29]  was  used. 
This  interpolation  is  based  on  the  linear  interpolation  along  the  edge  between  adjacent 
nodes.  The  shape  functions  are  rational  polynomials  that  interpolate  a  linear  held 
exactly.  For  completeness,  we  first  list  a  few  key  equations  from  [29],  and  then  derive 


90 


the  necessary  partial  derivatives  of  the  shape  functions. 

Let  n  be  the  number  of  vertices  of  the  polygon  and  i  be  an  index  of  those  n  sides. 
To  use  the  subscripts  to  indicate  the  pertinent  node,  the  coordinates  (xi,X2)  used  so 
far  will  be  renamed  (x,y)  in  this  section.  Each  edge  segment  l  is  described  by  a  line 
with  constant  coefficients  a  and  b. 


k  —  l  —  aa  -  b,;y 


(139) 


Each  polygon  will  have  the  center  node  as  its  local  origin.  The  coefficients  a  and  b 
are  then  determined  by  the  vertex  coordinates  in  the  reference  configuration. 
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The  line  coefficients  a  and  b  and  the  vertex  coordinates  are  used  to  calculate  the 
relative  weight  coefficients  k,  which  are  normalized  by  setting  k, i  to  a  value  of  one. 
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The  line  equations  and  coefficients  form  the  terms  of  -i/^,  the  numerator  of  the  shape 


function  associated  with  perimeter  node  i. 


j=n 

A(x,y)  =  Ki  JJ  lj(x,y)  (142) 

j^i+i 


The  denominator  polynomial  is  equal  to  the  sum  of  all  of  the  numerator  polynomi¬ 


als,  thus  forming  the  rational  polynomial  shape  function  Nt  for  each  node  i  on  the 
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perimeter  of  the  polygon. 


Ni(x,y)=4L-  (1«) 


Having  summarized  the  work  of  [29],  we  now  obtain  the  derivatives  of  the  shape 
functions  for  use  in  the  present  method.  Application  of  the  quotient  rule  leads  to  the 
expressions  for  the  shape  function  derivatives, 
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Each  of  the  terms  necessary  for  this  calculation  will  now  be  simplified.  Multiple 


applications  of  the  chain  rule  to  Eq.  (142)  result  in  the  expressions 
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Since  the  local  origin  was  placed  at  the  center  node,  and  the  center  node  is  the  only 
point  at  which  this  interpolation  will  be  applied,  all  of  the  non-constant  terms  of 
the  derivative  shape  functions  can  be  disregarded.  Then  the  calculations  for  the 
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individual  terms  of  Eqs.  (144)  and  (145)  reduce  to 
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These  formulas  are  substituted  into  Eqs.  (144)  and  (145)  to  provide  the  derivatives  of 
the  polygonal  shape  functions  at  the  center  nodes.  Thus  letting  Qt  be  scalar  values 
at  the  polygon  vertices  i  =  1  ...n,  the  partial  derivatives  of  Q  at  the  polygon’s  center 
node  are 
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~dy~  h 

To  improve  computational  efficiency,  much  like  global  stiffness  matrices,  the  indi¬ 
vidual  element  interpolation  functions  are  assembled  into  global  derivative  matrices 
Ta  =  d/dxa,  where  a  =  1,2  to  indicate  the  direction  of  the  partial  derivative.  To 
distinguish  between  the  sets  of  matrices,  let  the  superscripts  indicate  the  vector  upon 
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which  the  matrix  operates  (u  or  Q)  and  the  location  of  the  result  (c  for  triangle 
centroid  or  n  for  the  interior  node  of  a  polygon).  The  derivative  matrices  T"c  recover 
the  linear  triangle  partial  derivatives  from  the  nodal  displacements.  The  equations 
for  the  gradients  of  linear  triangles  may  be  found  in  most  introductory  finite  cle¬ 
ment  texts.  The  derivative  matrices  T“n  use  the  polygon  interpolations  to  recover 
the  partial  derivatives  of  u  at  the  polygon  interior  nodes  from  the  displacements  at 
the  perimeter  nodes.  Lastly.  Tj"  use  the  polygon  interpolations  to  recover  the  par¬ 
tial  derivatives  of  Q  at  the  polygon  interior  nodes  from  the  values  Q  at  the  triangle 
centroids. 

Residual  Calculations. 

The  residual  vector  R  contains  the  imbalance  between  the  internal  and  external 
forces  at  the  nodes  for  an  approximate  solution  of  the  governing  equations,  Eq.  (137). 
The  role  of  the  solver  is  to  reduce  the  size  of  the  residual  vector  to  an  acceptable 
level  with  the  user  supplying  the  metric  and  tolerance.  Using  the  entire  domain,  let 
ne  be  the  number  of  triangular  elements,  n  be  the  total  number  of  nodes,  and  nn  be 
the  number  of  interior  nodes.  The  following  quantities  are  calculated  sequentially  to 
return  the  residual  vector  R  to  the  solver. 

1.  Partial  derivatives  of  u  at  the  triangle  centroids  using  the  primal  mesh.  Six 
operations  are  required  (partials  of  three  u  vectors  in  two  directions)  using  the 
two  pre-calculated  discrete  derivative  operators  T"c. 

^  =  T“cu,  a  =  1,2,  *  =  1,2,3  (158) 

2.  Surface  unit  normal  vector  n  and  area  ratio  J a  at  the  center  nodes.  The  sur¬ 
face  normals  at  the  nodes  were  obtained  via  the  cross-product  of  two  tangent 
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vectors  [140,  15].  The  nodal  partial  derivatives  required  by  the  tangent  vectors 
can  be  obtained  by  gradient  recovery  techniques  or  by  applying  the  polygon 
interpolation  equations  (144)  and  (145)  to  the  full  element  patch.  The  latter 
technique  was  applied  using  the  two  pre-calculated  constant  matrices  Tun  in  six 
operations, 
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Thus  at  each  interior  node  the  tangent  vectors  were 
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Their  cross  product’s  magnitude,  JA,  is  the  ratio  of  the  deformed  area  to  the 
undeformed  area.  The  direction  provides  the  components  of  the  external  force 
vector.  The  unit  normal  vector  is 


gi  x  g2  _  gi  x  g2 
gi  X  gall  Ja 


(161) 


3.  External  force  f(u)  at  the  center  nodes.  Let  p  be  the  spatially-constant  mag¬ 
nitude  of  inflation  pressure  that  acts  normal  to  the  membrane  in  its  current 
configuration.  The  external  force  vector  is 


f  =  p  JAn 


(162) 


4.  Thickness  ratio  Jh  at  the  center  nodes.  Incompressibility  is  assumed  for  the 
deformed  thickness  calculation.  Since  the  volume  V  =  Ah  is  constant,  the 
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thickness  ratio  .//,  is  the  reciprocal  of  the  area  ratio. 


h_  _  V/A  _  Ao  _ 
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(163) 


5.  Qi  and  Q2  at  the  triangle  centroids  using  Eqs.  (135),  (136),  and  (138).  The 
deformed  thickness  h( u)  in  Eq.  (136)  will  need  to  be  replaced  by  Jhho-  Since 
Jh  will  be  calculated  at  the  center  node,  not  at  the  element  centroids  where  Q 
exists,  the  undeformed  membrane  thickness  h,0  is  used  for  now.  The  membrane 
thickness  will  be  corrected  in  the  last  step,  when  the  residual  vector  is  calculated. 

6.  Partial  derivatives  of  Qi  and  Q2  at  the  polygons’  center  nodes  by  applying  the 
two  pre-calcnlated  derivative  operators  T(-n  to  the  three  u  vectors. 

^=T«”Q^  a, /5  —  1,2  (164) 


7.  Global  residual  vector  R.  The  nodal  thickness  ratio  J/,  from  Eq.  (163)  is  now 
applied  to  the  partial  derivatives  of  Q  to  calculate  R,  G  9ft3,  residual  vector 
of  node  j.  The  nodal  residuals  are  assembled  into  the  global  residual  vector 
R  G  9ft3nn,  which  is  returned  to  the  nonlinear  solver  for  further  minimization. 
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In  the  point  collocation  method,  the  residual  equations  for  points  on  the  boundary 
are  determined  by  the  boundary  conditions  and  will  be  different  from  the  domain 
interior’s  governing  equation  [1],  Here,  for  the  homogeneous  boundary  conditions, 
the  displacements  are  simply  set  to  zero,  and  the  boundary  nodes  are  removed  from 
the  solution  vector. 
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Nonlinear  Solver. 


The  dynamic  simulation  under  development  which  will  use  this  membrane  model 
primarily  employs  an  accelerated  fixed-point  iteration  algorithm.  This  method  does 
not  require  the  construction  of  a  Jacobian  matrix  (or  even  an  approximation  of  one), 
so  significant  savings  in  computational  effort  are  expected.  Since  the  Newton- Raphson 
method  will  not  be  used,  the  tangent  stiffness  matrix  is  not  needed  and  will  not  be 
formulated  analytically.  Thus  for  the  static  cases  presented  here  it  was  necessary  to 
find  an  efficient  alternative  method  that  estimates  the  Jacobian  or  avoids  using  it 
altogether. 

For  this  study  we  utilized  “scipy .  optimize .  newton_krylov” ,  the  Newton- Krylov 
nonlinear  equation  solver  from  the  Python  library  SciPy  [65].  The  solver  was  supplied 
an  initial  guess  for  u  and  a  callable  function  that  returned  the  residual  vector  R(u), 
Eq.  (165).  The  loose  generalized  minimum  residual  (LGMRES)  method  was  selected 
as  the  inner  solver  [6]. 

Broadly  speaking,  Newton-Krylov  methods  employ  nested  iterative  solvers.  The 
outer  solver  performs  corrections  like  the  classical  Newton  method.  The  inner  solver 
is  one  of  many  linear  Krylov  subspace  methods  [131].  Jacobian-free  Newton-Krylov 
(JFNK)  methods  like  the  one  used  in  this  SciPy  routine  are  efficient  for  large  sys¬ 
tems  because  they  use  a  perturbation  of  the  entire  solution  vector  to  approximate 
Jacobian- vector  products;  this  approach  is  more  efficient  than  the  finite  difference 
Jacobian  construction,  which  requires  a  perturbation  of  each  element  of  the  solution 
vector  [75,  69].  The  actual  Jacobian  is  not  needed,  yet  convergence  can  approach  that 
of  the  Newton  method.  Preconditioning  by  providing  an  approximation  to  the  Jaco¬ 
bian  is  highly  recommended  and  often  necessary  for  adequate  performance  [12,  131]. 
However,  for  the  simulations  in  this  study,  performance  was  more  than  adequate 
without  supplying  a  preconditioner  to  the  solver. 


Because  stress  stiffening  is  the  source  of  transverse  resistance,  the  solver  will  fail 
when  starting  from  a  flat,  slack  membrane.  There  are  several  remedies  for  this  prob¬ 
lem.  The  approach  used  in  this  study  was  to  prescribe  an  initially  non-flat  shape 
([141]  used  this  approach  for  a  box-shaped  membrane  and  mentioned  its  necessity). 
Simple  parabolic  profiles  were  sufficient.  A  second  option,  dynamic  relaxation,  uses  a 
dynamic  model  with  damping  to  settle  to  the  static  solution  [64,  28,  140].  This  option 
is  convenient  since  it  is  not  necessary  to  code  a  separate  solver  for  the  static  solution; 
and  it  also  has  the  benefit  of  verifying  some  aspects  of  the  dynamics  code.  How¬ 
ever,  convergence  can  be  extremely  slow  if  the  damping  mechanism  is  not  carefully 
designed.  Some  of  these  challenges  can  be  avoided  by  using  pseudo-transient  continu¬ 
ation  [40,  70].  The  non-physical  time  step  sizes  can  be  controlled  using  the  successive 
evolution-reaction  (SER)  technique  [75].  A  third  approach,  also  not  employed  herein, 
entails  initially  applying  a  pretension  to  aid  the  solver  and  later  removing  it  for  the 
final  solution. 

It  is  worth  mentioning  a  few  other  numerical  alternatives  that  may  be  used  in  the 
absence  of  an  analytical  Jacobian.  A  more  comprehensive  review  of  these  alternatives 
may  be  found  in  [140].  Finite  difference  Jacobian  approximations  are  easy  to  perform 
but  are  not  robust,  and  they  become  very  expensive  as  the  size  of  the  problem  in¬ 
creases.  Automatic  differentiation  extracts  the  derivatives  directly  from  the  code  [47] . 
Multigrid  methods  [21,  52]  (in  particular  the  Full  Approximation  Scheme  [20]),  and 
combinations  of  Newton-Krylov  and  multigrid  methods  [66,  41]  require  some  careful 
coding,  but  the  ultimate  computational  efficiency  gains  can  be  impressive.  Lastly, 
depending  upon  how  the  method  is  formulated,  dynamic  relaxation  may  also  be  a 
viable  alternative. 
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Verification. 


Patch  tests  were  performed  for  the  linear  plane  stress  scenario  [58].  All  of  the  tests 
featured  a  unit  vertical  rigid  body  translation.  The  square  patch  measured  two  units 
per  side  and  was  centered  at  (1, 1)  as  shown  in  Figure  29.  The  patch  contained  a  single 
interior  node  at  (1.6, 1.4),  approximately  on  the  perimeter  of  its  dual-mesh  polygon. 
Given  a  linear  displacement  field,  boundary  nodes  were  displaced  accordingly  and 
the  displacement  of  the  center  node  was  checked  against  the  exact  field.  Introducing 
Young’s  modulus  E,  the  physical  constants  were  set  to  E  =  1000,  h  —  1,  and  v  =  0.5. 
The  results  in  Table  7  show  that  the  model  exactly  reproduced  the  constant  strains 
and  stresses,  even  with  irregular  elements. 

Further  verification  and  convergence  determination  was  performed  using  the  Method 


Figure  29.  Square  2x2  5-node  patch  for  the  patch  test.  Solid  lines  are  the  triangles  of 
the  primary  mesh;  dashed  lines  illustrate  the  node-centered  polygon  of  the  dual  mesh. 

of  Manufactured  Solutions  [112,  23].  In  this  method,  a  solution  is  fabricated  (it  does 
not  have  to  be  physically  plausible),  and  the  forcing  function  is  calculated  from  the 
governing  PDE.  Then,  the  forcing  function  is  used  in  the  numerical  model  to  obtain  an 
approximate  solution.  The  error  is  then  the  difference  between  the  model’s  solution 
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Table  7.  Patch  test  results.  A  star  (*)  indicates  a  non-zero  magnitude  of  less  than 
10~n;  strains  and  stresses  for  all  four  elements  were  reported  as  exact  when  integer 
values,  and  exact  to  at  least  eight  decimal  places  when  fractional. 


u  held 

u(1.6,1.4) 

E 

N jh 

(0,1,1) 

(1,0*, 1) 
(0*,1,1) 

(0,0*, 0*) 
(0,0* ,0*) 

*  ~  ~ 

o  o 

t-H  t-H 

cT  j*T 

o. 

(1.6,0*, 1) 
(0*, 1.6,1) 

(1.5,0*,0) 
(0.5,0*, 0.5) 

(2000, 1000,  0) 
(666.7,333.3,333.3) 

(y,o,i) 

(o,y,i) 

(1.4,0*, 1) 
(0*, 1.4,1) 

(0,0.5, 0.5) 
(0, 1.5,0*) 

(333.3,666.7,333.3) 
(1000,  2000,  0*) 

and  the  manufactured  solution.  The  manufactured  solution,  Eq.  (166),  was  devised 
in  accordance  with  the  recommended  guidelines  found  in  [115].  Non-unity  constants 
were  chosen  such  that  the  solution  magnitudes  and  derivatives  were  of  approximately 
the  same  order  of  magnitude  in  all  three  axes  such  that  potential  formula  errors  might 
be  revealed. 


u  (xi,  x2)  =  { 


0.1  Te^1/2-12/2) 
-0.37e(-Xl/4~X2/2) 
0.71e(Xl+Xa) 


(166) 


The  lengthy  formulas  for  the  corresponding  anisotropic  pressure  vector  were  obtained 
by  substituting  the  manufactured  solution  u  of  Eq.  (166)  into  Eqs.  (132),  (133),  (134), 
(135),  and  (136).  This  process  was  performed  independently  of  the  point  collocation 
code.  The  pressure  vector  formulas  were  generated  symbolically  with  the  Python 
library  SyrnPy  [129]  and  inserted  into  the  point  collocation  model.  The  Dirichlet 
boundary  conditions  were  satisfied  by  constraining  perimeter  node  displacements  in 
accordance  with  the  manufactured  solution,  Eq.  (166). 

To  investigate  the  effects  of  discretization  error,  verification  was  performed  on 
two  domains:  a  hexagon  and  a  circle.  The  hexagonal  domain  was  ideally  discretized 
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with  a  structured,  symmetrical  primary  mesh  of  equilateral  triangles  to  minimize 
discretization  error.  The  six  points  of  the  domain  laid  on  the  unit  circle.  The  unit- 
radius  circular  domain  was  discretized  with  asymmetrical,  unstructured  grids — the 
same  grids  which  will  be  used  later  in  the  validation  phase.  Three  grids  from  each 
domain  are  shown  in  Figure  30.  All  grids  were  created  using  the  open  source  software 
Gmsh  [43],  which  contains  a  “refine  by  splitting”  feature  (elsewhere  called  “refine  by 
quartering”  [74])  to  easily  perform  structured  refinement  with  a  grid  ratio  of  two. 
For  this  study  the  representative  element  size  of  a  mesh,  he,  was  calculated  as  the 
length  of  the  side  of  an  equilateral  triangle,  where  the  area  of  the  triangle  was  equal 
to  the  mean  element  area.  The  error  measures  eu  and  6e  were  the  Euclidean  norms  of 
the  displacement  error  vector  and  the  strain  error  vector  at  the  origin,  respectively. 
The  strain  was  recovered  as  the  mean  of  the  strains  of  the  surrounding  triangular 
elements.  Both  measures  were  normalized  by  the  value  of  the  exact  solution. 


Figure  30.  Point  collocation  static  verification  and  validation  meshes.  Top  row:  Hexag¬ 
onal  grids  with  he  =  0.5,  0.25,  and  0.0625  (hex2,  hex4,  and  hexl6).  Bottom  row:  The 
three  finest  circular  grids,  used  for  both  verification  and  validation  (circle3,  circle2,  and 
circlel,  left-to-right). 


The  convergence  study  results  are  shown  in  Table  8  and  Figure  31.  The  observed 


102 


order  of  accuracy  p  converged  to  two  for  both  displacement  and  strain  in  the  hexag¬ 
onal  domain.  The  displacement  errors  were  slightly  higher  for  the  circular  grids  due 
to  the  discretization  error  of  the  unstructured  grid;  however,  the  same  order  of  con¬ 
vergence  was  observed.  The  strain  convergence  in  the  circular  grid  did  not  behave 
as  neatly  because  the  individual  strain  vector  components  converged  differently:  Eu 
from  above,  EV2  from  below,  and  E22  non-monotonically.  Thus  the  vector  norm  con¬ 
verged  non-monotonically.  The  slope  of  a  least-squares  linear  fit  of  the  four  circular 
grids’  strain  errors  provided  p  =  2.24,  more  in  line  with  the  other  convergence  rates. 
Interestingly,  the  circular  grids  produced  more  accurate  strain  predictions  than  the 
hexagonal  grids  despite  the  non-monotonic  strain  convergence  behavior  and  lower 
displacement  accuracy. 

Table  8.  Convergence  study  results.  The  hexagonal  grids  are  labeled  with  the  number 
of  elements  from  the  origin  to  the  vertex  at  (1,0).  The  circular  grids  are  numbered 
sequentially  from  fine  to  coarse. 


Grid 

he 

log  (he) 

log(e„) 

Pu 

log  (eg) 

Pe 

hex2 

0.5 

-0.30 

-1.04 

- 

-0.98 

- 

hex4 

0.25 

-0.60 

-1.74 

2.32 

-1.75 

2.56 

hex8 

0.125 

-0.90 

-2.36 

2.07 

-2.39 

2.13 

hexl6 

0.0625 

-1.20 

-2.97 

2.02 

-3.00 

2.03 

hex32 

0.03125 

-1.51 

-3.57 

2.00 

-3.60 

2.01 

hex64 

0.015625 

-1.81 

-4.17 

2.00 

-4.21 

2.00 

circle4 

0.537 

-0.27 

-0.77 

- 

-1.02 

- 

circle3 

0.273 

-0.56 

-1.47 

2.39 

-2.93 

6.50 

circle2 

0.137 

-0.86 

-2.10 

2.10 

-2.80 

-0.42 

circlel 

0.069 

-1.16 

-2.72 

2.05 

-3.29 

1.63 

The  verification  process  has  demonstrated  that  the  model  correctly  and  consis¬ 
tently  solved  the  coded  governing  equations.  The  method’s  observed  order  of  conver¬ 
gence  of  both  displacement  and  strain  was  two.  Also,  the  circular  grids  introduced 
an  acceptably  small  amount  of  discretization  error  and  are  therefore  suitable  for  use 
in  the  validation  phase. 
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Figure  31.  Convergence  study.  The  triangle  displays  a  reference  slope  of  two.  Dashed 
lines  are  hexagonal  grid  results;  solid  lines  are  circular  mesh  results.  Grids  hex2  and 
circle4  from  Table  8  are  omitted. 


Validation. 

The  point  collocation  model  was  validated  against  the  experimental  and  finite 
element  model  results  in  [125].  In  that  study,  a  latex  rubber  sheet  was  inflated  from 
below,  and  the  displacement  field  was  extracted  by  optically  tracking  a  random  speck¬ 
ling  pattern  on  the  membrane  surface.  Strains  were  calculated  from  the  displacement 
field  during  post-processing.  A  geometrically  nonlinear  finite  element  model  was  com¬ 
pared  to  another  finite  element  approximation  by  [122],  then  validated  against  the 
experimental  data. 

Two  cases  were  selected  for  validation  in  this  paper:  one  with  and  the  other 
without  uniform  prestress.  In  the  case  without  prestress,  model  predictions  were 
compared  directly  to  the  experimental  data.  The  physical  constants  were  outer  radius 
R  =  57.15  mm,  thickness  h  =  0.12  mm,  modulus  of  elasticity  E  =  2  MPa,  and  Poisson 
ratio  v  =  0.5. 

In  the  prestressed  case,  a  non-isotropic  prestrain  field  precluded  use  of  the  ex¬ 
perimental  results.  The  geometrically  nonlinear  finite  element  model  was  therefore 
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used  as  a  benchmark  for  this  study.  The  FEM  model  was  suitable  as  a  benchmark 
because  it  was  validated  in  [125]  against  the  experimental  data  for  the  same  problem 
configuration  (including  geometry,  boundary  conditions,  material,  and  load  type)  and 
to  much  greater  loads  and  deformations.  The  physical  constants  for  this  case  were 
R  =  3.5  mm,  h  =  1.0/mi,  E  =  71.0  GPa,  and  v  =  0.345.  Following  convention  for  the 
Hencky  problem,  the  lateral  deflection  is  normalized  as  w/R,  and  the  nondimensional 
pressure  q  normalizes  the  inflation  pressure  as  q  =  pR/ Eh. 

Representative  grid  convergence  results  are  shown  in  Table  9.  The  format  of  the 
presentation  comes  from  [24] ,  which  put  forth  useful  guidelines  to  standardize  the  re¬ 
porting  of  CFD  numerical  study  results.  As  explained  earlier,  all  grid  ratios  are  equal 
to  two.  The  symbol  <p  represents  the  magnitude  of  the  field  variable  (displacement  or 
strain).  A  subscript  1,  2,  or  3  indicates  the  results  were  obtained  from  the  grid  circlet, 
circle2,  or  circle3  respectively  (see  Table  8  and  Figure  31).  The  calculated  order  of 
accuracy  is  p.  The  symbol  indicates  the  value  was  obtained  from  Richardson 
extrapolation  using  grids  2  and  3.  The  approximate  and  extrapolated  relative  errors 
(magnitude  percentage  change  in  (j)  from  one  grid  to  the  next  finer)  are  shown  as  ea 
and  text' 

The  convergence  study  confirms  that  the  grids  were  sufficiently  refined,  as  indi¬ 
cated  by  the  convergence  of  the  extrapolated  values  <j)ext  and  the  small  relative  errors 
e.  The  extrapolated  values  were  taken  as  the  model  solution  for  the  remain¬ 
ing  discussions.  The  second-order  observed  accuracy  from  the  verification  phase  was 
maintained  for  both  displacement  and  strain  within  the  range  of  the  model’s  expected 
applicability.  The  order  of  convergence  of  w/R  deteriorated  to  below  one  for  large 
deformations  (i.e.,  for  w/R  >  0.4). 

The  point  collocation  model  predictions  are  compared  to  the  experimental  results 
in  Figure  32  and  to  the  FEM  model  benchmark  in  Figure  33.  Agreement  is  excellent 
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Table  9.  Validation  convergence  study  for  displacement  ( w/R )  and  radial  strain  ( Er )  at 
the  origin. 


Dep.  variable  0 

w/R 

w/R 

Er 

Er 

w/R 

Prestress 

None 

None 

None 

None 

250MPa 

q 

0.024 

0.048 

0.048 

0.083 

0.0030 

03  (coarse  mesh) 

0.1748 

0.2223 

0.0339 

0.0525 

0.08005 

02  (medium) 

0.1739 

0.2212 

0.0348 

0.0541 

0.07983 

0i  (fine) 

0.1737 

0.2210 

0.0350 

0.0546 

0.07978 

P 

2.00 

1.96 

2.04 

1.98 

2.05 

cb32 

t* ext 

0.1736 

0.2209 

0.0350 

0.0547 

0.07976 

621 
r  ext 

0.1736 

0.2209 

0.0350 

0.0547 

0.07976 

e21 

0.13% 

0.12% 

0.59% 

0.76% 

0.07% 

e21 

ext 

0.04% 

0.04% 

0.19% 

0.26% 

0.02% 

with  and  without  prestress  up  to  w/R  ~  0.25,  beyond  which  the  point  collocation 
model  begins  to  underestimate  the  displacement.  This  behavior  is  in  agreement  with 
the  comparison  of  Mooney- Rivlin  and  Hookean  material  models  in  [106],  so  the  grad¬ 
ual  loss  of  accuracy  beyond  this  point  may  be  attributed  to  the  onset  of  hyperelastic 
effects.  The  nodal  strain  predictions  continue  to  match  the  experimental  data  beyond 
w/R  ~  0.45.  Just  as  in  the  verification  phase,  the  strains  are  actually  more  accurate 
than  the  displacements.  Typically  the  directly-calculated  derived  variables  (strains 
and  stresses)  converge  more  slowly  than  the  displacements,  though  post-processing 
recovery  procedures  may  improve  the  accuracy  [146,  98].  The  staggered  mesh  of  this 
model  uses  constant-strain  triangles,  so  nodal  strain  values  must  be  recovered  by  one 
of  the  gradient  recovery  procedures.  Simple  averaging  from  neighboring  elements  was 
sufficient  to  produce  the  excellent  relative  accuracy  of  strain  at  the  center  node. 

Validation  against  the  experimental  data  has  shown  that  the  model  is  accurate 
for  displacements  up  to  w/R  ~  0.25,  at  which  point  a  hyperelastic  model  would  be 
more  appropriate  as  discussed  in  [125]. 
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Figure  32.  Validation  with  no  prestress  by  comparison  to  experimental  data. 
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Figure  33.  Validation  of  the  prestressed  case  by  comparison  to  the  FEM  benchmark. 
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Summary. 


A  membrane  model  intended  for  eventual  use  in  dynamic  aeroelastic  simulations 
was  presented  in  this  section,  and  its  performance  for  membranes  at  static  equilib¬ 
rium  was  investigated.  The  model  effectively  combines  several  unconventional  for¬ 
mulations  in  structural  engineering,  including  a  staggered  grid  with  robust  low-order 
interpolation  schemes,  grouped  nonlinear  products  as  degrees  of  freedom,  and  the 
point  collocation  method.  Method  capabilities  include  variable  thickness,  follower 
forces,  and  arbitrary  prestress.  Rigorous  verification  demonstrated  consistency,  and 
the  observed  order  of  convergence  was  two  for  both  displacement  and  strain.  Dur¬ 
ing  validation  with  respect  to  a  static  circular  membrane  (the  Hencky  problem),  the 
point  collocation  model  predictions  agreed  with  experimental  data  and  benchmark 
FEM  code  until  the  region  where  hyperelastic  response  began  to  dominate. 

The  primary  feature  that  distinguishes  this  approach  is  its  simplicity.  Element 
integration  is  avoided  entirely.  The  group  formulation  permits  the  same  treatment 
in  all  three  axes,  and  the  resulting  code  is  explicit  and  self-documenting.  Overall, 
the  framework  of  the  approach  is  highly  modular  and  flexible.  Any  given  step  can  be 
performed  by  interchangeable  subroutines.  For  example,  the  polygon  interpolation 
technique  could  be  replaced  by  least  squares  or  radial  basis  function  routines  without 
replacing  the  remaining  code.  The  residual  subroutine  readily  accepts  different  strain- 
displacement  or  material  models  (including  nonlinear  models). 

Membrane  Dynamics 

This  section  presents  the  model  that  fulfills  the  objectives  of  this  research  by 
accurately  simulating  the  transient  response  of  a  geometrically  nonlinear  membrane. 
The  model  incorporates  the  new  temporal  and  spatial  discretization  schemes,  namely 
the  Hermite  time  interpolation  method  of  Section  IV  and  the  point  collocation  method 
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of  Section  IV. 


The  complete  dynamic  membrane  model  will  be  verified  using  the  Method  of 
Manufactured  Solutions.  Rates  of  convergence  will  be  calculated  from  the  results. 
A  free-response  case  study  will  be  performed  using  the  converged  results  of  the  two 
time-integration  methods  as  the  baseline  solution.  Throughout  this  section,  results 
of  the  model  using  the  Hermite  integration  method  will  be  compared  to  results  using 
the  second-order  Newmark  trapezoid  method  (shortened  to  “the  Newmark  method” 
from  now  on).  Despite  the  disparity  in  the  orders  of  the  methods,  the  Newmark 
method  was  chosen  due  to  its  prevalence  and  familiarity  in  the  structural  engineering 
community.  Additionally,  since  two  different  dynamic  time  scales  were  observed,  a 
hybrid  method  will  be  investigated  that  uses  the  Newmark  trapezoid  method  for 
in-plane  dynamics  and  the  Hermite  interpolation  method  for  out-of-plane  dynamics. 

The  governing  equations  for  the  dynamic  nonlinear  membrane  were  derived  in  [111]. 
They  are  identical  to  those  of  Section  IV  (Eq.(132),  Eq.(133),  and  Eq.(134))  but  with 
the  inertial  forces  retained.  The  density  p  is  in  terms  of  mass  per  unit  volume. 


phiii  = 
phu2  = 
phu3  = 


d  Nu  dN12 
dxi  dx2 
dN22  |  8N12 
dx2  dxi 

JL(n  ^ 

dx\  \  11  dx\ 


+  /i 


+  /2 


+  A 12 


d 

dx2 


du3 

dx2 


+  N12 


+  h 


(167) 

(168) 
(169) 


Retaining  from  Section  IV  the  expressions  for  strains  and  stresses,  as  well  as 
the  definitions  for  the  vectors  Q,  the  vectorized  version  of  the  dynamic  governing 
equations  is  therefore 


phii 


<9Qi  dQ2 

dxi  dx2 


(170) 
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Jerk  Constraints. 


The  jerk  constraints  were  derived  by  taking  the  time  derivative  of  the  equilibrium 
equations,  Eq.  170,  to  obtain 


j 


J_  (  .  dQa  • 

ph  1  dx\  dx2 


The  jerk  constraints  require  formulas  for  the  strain  rates  of  change, 


(171) 


En 

E22 

2Ei2 


dv\  f  du\  dvi 


dxi 

dv2 

dx2 

dv2 


+ 


+ 


dxi  dx\ 
dui  dvx 


dxo  dx- 


+ 


+ 


du2  dv2 
dxi  dxi 
du2  dv2 
dxo  dxo 


+ 


+ 


du3  dv3 
dxi  dxi 
du3  dv3 
dxo  dxo 


dv\ 
dxi  dx2 
dv2  du2 


+ 


dxi  dx2 


+ 


dv\  du\ 
dx\  dx2 
du2  dv2 
dxi  dx2 


+ 


dui  dv\ 


+ 


dxi  dx2 
dv3  du3 


dx\  dx2 


+ 


du3  dv3  \ 

dx\  dx2  J 


(172) 

(173) 

(174) 

(175) 


which  after  following  the  same  process  as  described  above  provide  the  rates  of  change 

of  Q, 


Qi 


Nn 

N12 


V 


Q2 


N12 

N22 


> 


(176) 


(177) 


The  calculation  of  the  jerk,  j,  uses  much  of  the  same  code  as  that  of  the  acceleration,  a. 
The  spatial  discretization  algorithms  and  data  structures  are  utilized  for  both,  as  well 
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as  surface  gradients.  Thus,  in  terms  of  computational  expense,  where  the  calculation 
of  a  is  one  function  call,  the  work  required  to  calculate  both  a  and  j  is  less  than  two 
function  calls.  As  detailed  in  Section  IV,  from  classical  analysis  of  a  linear  oscillator, 
the  resulting  scheme  is  conditionally  stable,  frequency-preserving,  non-dissipative, 
and  more  accurate  than  the  Newmark  trapezoidal  and  central  difference  methods  by 
several  orders  of  magnitude. 

Verification. 

Performance  analysis  was  accomplished  using  the  Method  of  Manufactured  So¬ 
lutions  [112,  23].  In  this  procedure,  a  solution  is  fabricated  (it  does  not  have  to  be 
physically  plausible)  and  the  forcing  function  is  analytically  derived  from  the  govern¬ 
ing  PDE  and  the  constitutive  equations.  The  forcing  functions,  having  been  derived 
through  symbolic  differentiation,  are  discretization-independent.  When  the  simula¬ 
tions  are  run  for  a  discretized  domain,  the  Dirichlet  boundary  conditions  are  set  by  the 
manufactured  solution  and  the  derived  forcing  function  is  applied  across  the  domain. 
The  error  is  then  the  difference  between  the  model’s  solution  and  the  manufactured 
solution. 

The  domain  was  defined  to  be  a  hexagon  circumscribed  by  a  circle  of  radius 
R  =  57.15  mm  (the  radius  of  the  circular  membrane  in  Ref.  [125]).  See  Fig.  34  for  a 
depiction.  The  origin  was  placed  at  the  center  of  the  domain.  This  domain  contained 
the  structured  primary  mesh,  which  consisted  of  equilateral  triangles  with  side  length 
he  =  R/8. 

The  manufactured  solution  crafted  for  this  study,  Eq.  (178),  dictated  a  smooth, 
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Figure  34.  Hexagonal  grid  for  the  dynamic  membrane  verification,  showing  the  cir¬ 
cumscribed  circle  of  radius  R. 


curved  surface  by  the  displacement  fields, 


u(x,2/,t)  =  A(t)  < 


3y/2 


3x+y 


(178) 


where  the  time-varying  amplitudes  were 


A(t) 


R 

To 


(cos500f  +  2) 


(179) 


At  the  origin,  the  exact  solutions  for  displacement  in  all  three  axes  were  precisely 
A(t),  so  it  was  known  a  priori  that  any  frequency  content  where  u  ^  500  rad/sec 
was  induced  by  model  errors.  The  physical  constants  were  chosen  to  be  those  of  the 
experimental  membrane  in  Ref.  [125]:  membrane  thickness  h  =  0.12  mm,  modulus  of 
elasticity  E  =  2  MPa,  density  p  =  990  kg/m3,  and  Poisson  ratio  v  =  0.5.  The  lengthy 
formulas  for  the  pressure  vectors  and  their  rates  of  change  that  produce  the  solution 
of  Eq.  166  were  generated  symbolically  with  the  Python  library  SymPy  [129] . 

A  convergence  study  was  performed  to  compare  the  Hermite  interpolation  method 
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with  physical  jerk  constraints  from  Eq.  171  to  the  Newmark  trapezoid  method,  where 
7  =  0.5  and  f3  =  0.25.  The  simulation  was  run  for  one  period  of  the  exact  solution 
( T  =  27t/500  s)  using  20,  40,  80,  and  160  time  steps.  First,  the  results  with  respect  to 
the  final  displacement  errors  at  the  origin  will  be  discussed,  followed  by  results  with 
respect  to  mean  errors  across  all  of  the  spatial  nodes  throughout  the  simulation. 

The  convergence  data  using  the  final  displacements  u i ,  u2,  and  m3  at  the  origin  are 
listed  in  Table  10  and  plotted  in  Fig.  35,  with  the  out-of-plane  errors  shown  separately 
from  the  in-plane  errors.  Error  was  normalized  as  a  percentage  of  the  exact  solution’s 
amplitude.  The  values  shown  are  log10  | error |.  An  asterisk  indicates  invalid  data  due 
to  non-monotonic  convergence,  and  “n/a”  indicates  lack  of  data  due  to  failure  of  the 
algorithm  to  converge.  The  rates  of  convergence,  p,  are  shown  for  each  increment 
of  mesh  refinement.  The  uneven  trend  lines  of  the  in-plane  errors  for  the  trapezoid 
method  were  produced  by  their  non-monotonic  convergence.  The  Hermite  method 
improved  upon  the  Newmark’s  accuracy  by  three  orders  of  magnitude  for  out-of-plane 
errors,  but  the  techniques  were  comparable  for  in-plane  displacements.  The  observed 
rate  of  convergence  was  approximately  two  for  both  methods. 

Table  10.  Dynamic  nonlinear  membrane  verification  convergence  data  for  the  final 
displacements  at  the  origin. 


Trapezoid 

Hermite 

Steps 

Ml 

u2 

u3 

Ml 

u2 

u3 

20 

-0.60 

-0.61 

0.98 

n/a 

n/a 

n/a 

40 

-2.38 

-1.93 

0.36 

-1.98 

-1.98 

-2.72 

80 

-1.86 

-1.95 

-0.25 

-2.52 

-2.52 

-3.31 

160 

-2.39 

-2.46 

-0.86 

-3.06 

-3.05 

-4.06 

p  20  — y  40 

5.90* 

4.38* 

2.06 

n/a 

n/a 

n/a 

p  40  — >  80 

-1.71* 

0.07* 

2.03 

1.79 

1.78 

1.98 

p  80  — » 160 

1.74 

1.70 

2.01 

1.80 

1.78 

2.48 

The  convergence  study  results  using  the  mean  error  magnitude  for  all  spatial 
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Time  steps  per  period 


(a)  Out  of  plane  (b)  In-plane 

Figure  35.  Dynamic  nonlinear  membrane  convergence  plots  for  the  final  displacements 
at  the  origin. 


nodes  over  all  time  steps  are  shown  in  Table  11  and  plotted  in  Fig.  36.  The  error  is 
defined  as  error  =  ^  Uj  —  Uj\/Nn(N  +  1),  where  u  are  the  exact  solutions,  Nn  is  the 
number  of  nodes  in  the  mesh,  and  N  is  the  number  of  time  steps  in  the  simulation. 
The  values  shown  are  log10  |  error  j.  The  entry  “n/a”  indicates  lack  of  data  due  to 
failure  of  the  algorithm  to  converge.  The  rates  of  convergence,  p,  are  shown  for  each 
increment  of  mesh  refinement.  The  observed  rate  of  convergence  for  both  methods 
was  two,  the  same  as  that  of  the  final  displacement  errors  at  the  origin  in  Fig.  35(a). 
With  this  metric,  the  accuracy  of  the  Hermite  method  was  approximately  1.5  orders 
of  magnitude  better  than  the  trapezoid  method. 

The  reason  why  the  Hermite  method’s  accuracy  was  significantly  greater  for  out- 
of-plane  displacements  than  for  in-plane  displacements  is  depicted  clearly  in  Fig.  37, 
a  representative  illustration  of  the  errors  over  time  from  a  separate  simulation  using 
the  trapezoid  method.  The  representative  plot  of  normalized  model  errors  was  gener¬ 
ated  using  the  Newmark  trapezoid  method  and  external  forcing  for  a  manufactured 
solution.  It  is  evident  that  the  oscillation  frequency  of  in-plane  errors  significantly  ex- 
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Table  11.  Dynamic  nonlinear  membrane  convergence  data  for  the  mean  error  magni¬ 
tudes  throughout  the  simulation. 


Steps 

Trapezoid 

Hermite 

20 

0.25 

n/a 

40 

-0.35 

-1.98 

80 

-0.96 

-2.57 

160 

-1.56 

-3.17 

p  20  — >  40 

2.00 

n/a 

p  40  — »  80 

2.00 

1.98 

p  80  — >  160 

2.00 

1.98 

Figure  36.  Dynamic  nonlinear  membrane  convergence  plot  of  the  mean  error  magni¬ 
tudes  for  all  three  axes,  all  spatial  nodes,  all  time  steps. 
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ceeded  that  of  the  out-of-plane  error.  A  membrane  is  much  stiffer  in-plane,  where  the 
tensile  internal  forces  respond  linearly  with  displacement.  In  contrast,  the  resistance 
to  out-of-plane  motion  is  zero  for  a  flat  membrane,  and  stress-stiffening  occurs  with 
ever  greater  transverse  deformation.  Thus  the  present  system  has  two  time  scales  for 
an  integrator  to  deal  with.  The  Hcrmite  method’s  accuracy  can  be  fully  realized  for 
the  slow  signal,  where  it  significantly  outperforms  the  trapezoid  method.  However, 
the  fast  in-plane  response  can  produce  oscillatory  results  or  cause  the  conditionally 
stable  Hermite  method  to  diverge  prior  to  run  completion  as  the  time  step  increases, 
and  in  fact  the  solution  did  not  converge  for  the  time  step  size  A t/T  =  0.05.  (Note 
that  the  period  T  in  this  calculation  is  the  period  of  the  manufactured  solution,  not 
of  the  numerical  solution’s  oscillations  as  depicted  in  Fig.  37.) 


Figure  37.  Illustration  of  the  nonlinear  membrane’s  two  time  scales:  rapid  in-plane 
responses  {u-\  and  M2),  and  slower  out-of-plane  dynamics  (M3). 


Following  these  observations,  a  hybrid  technique  was  explored  to  capitalize  on 
the  strengths  of  both  the  Newmark  and  the  Hermite  methods:  the  unconditional 
stability  of  the  trapezoid  method  for  the  rapidly  changing  in-plane  displacements, 
and  the  significantly  improved  accuracy  of  the  Hermite  method  for  the  slower  out-of¬ 
plane  dynamics.  It  was  proved  in  Section  IV  that  the  Newmark  could  be  obtained 
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by  setting  the  jerk  constraints  appropriately.  Therefore,  coding  was  straight-forward 
to  utilize  the  different  jerk  constraints  for  in-plane  versus  out-of-plane  equations  of 
motion. 

The  results  of  the  previous  convergence  studies  are  now  presented  again,  but  with 
the  hybrid  method  included  for  comparison.  The  convergence  data  using  the  final 
displacements  at  the  origin  are  listed  in  Table  12  and  plotted  in  Fig.  38  (with  data 
carried  forward  from  Table  10  and  Fig.  35).  The  convergence  study  results  using  the 
mean  error  magnitude  for  all  spatial  nodes  over  all  time  steps  are  shown  in  Table  13 
and  plotted  in  Fig.  39  (with  data  carried  forward  from  Table  11  and  Fig.  36). 

As  expected,  the  convergence  data  show  that  the  hybrid  method  effectively  com¬ 
bined  the  strengths  of  the  Hcrmite  and  Newmark  methods.  In-plane  errors  were 
comparable  to  the  pure  trapezoid  method,  but  the  out-of-plane  accuracy  was  bet¬ 
ter  than  the  trapezoid  method  by  two  orders  of  magnitude  as  shown  in  Fig.  38(a). 
Using  the  mean  error  norm,  the  hybrid  method  errors  were  an  order  of  magnitude 
smaller  than  those  of  the  trapezoid  method.  In  addition,  as  far  as  the  improving  sta¬ 
bility,  the  hybrid  method  successfully  completed  the  simulation  with  time  step  size 
A t/T  =  0.05,  where  the  purely  Hermite  method  failed. 

Table  12.  Dynamic  nonlinear  membrane  verification  convergence  data  for  the  final 
displacements  at  the  origin,  including  the  hybrid  method. 


Trapezoid 

Hermite 

Hybrid 

Time  Steps 

Ml 

U-2 

u3 

Ml 

u2 

U3 

Mi 

u2 

u3 

20 

-0.60 

-0.61 

0.98 

n/a 

n/a 

n/a 

-0.55 

-0.56 

-0.88 

40 

-2.38 

-1.93 

0.36 

-1.98 

-1.98 

-2.72 

-1.06 

-1.04 

-1.53 

80 

-1.86 

-1.95 

-0.25 

-2.52 

-2.52 

-3.31 

-1.90 

-1.86 

-2.18 

160 

-2.39 

-2.46 

-0.86 

-3.06 

-3.05 

-4.06 

-2.55 

-2.51 

-2.80 

p  20  — >  40 

5.90* 

4.38* 

2.06 

n/a 

n/a 

n/a 

1.68 

1.60 

2.16 

p  40  — *  80 

-1.71* 

0.07* 

2.03 

1.79 

1.78 

1.98 

2.79 

2.74 

2.15 

p  80  — »  160 

1.74 

1.70 

2.01 

1.80 

1.78 

2.48 

2.17 

2.16 

2.06 
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Figure  38.  Dynamic  nonlinear  membrane  convergence  plots  for  the  final  displacements 
at  the  origin,  including  the  hybrid  method. 


Table  13.  Dynamic  nonlinear  membrane  convergence  data  for  the  mean  error  magni¬ 
tudes  throughout  the  simulation,  including  the  hybrid  method. 


Time  Steps 

Trapezoid 

Hermite 

Hybrid 

20 

0.25 

n/a 

-0.84 

40 

-0.35 

-1.98 

-1.43 

80 

-0.96 

-2.57 

-2.04 

160 

-1.56 

-3.17 

-2.65 

p  20  — >  40 

2.00 

n/a 

1.97 

p  40  — »  80 

2.00 

1.98 

2.01 

p  80  — »  160 

2.00 

1.98 

2.05 
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Figure  39.  Dynamic  nonlinear  membrane  convergence  plot  of  the  mean  error  magni¬ 
tudes  for  all  three  axes,  all  spatial  nodes,  all  time  steps,  including  the  hybrid  method. 

Free  Response  Analysis. 

A  case  study  of  the  undamped  free  response  of  a  membrane  was  performed  to 
challenge  the  model  with  a  rapid  initial  transient  and  wide  range  of  response  frequen¬ 
cies.  The  initial  at-rest  configuration  was  a  smooth  peak  slightly  off  the  center  of  a 
square  membrane  (see  Eq.  180  and  Fig.  40(a)),  designed  to  introduce  traveling  waves. 
Initial  in-plane  displacements  were  zero. 

u3(Xl,  x2)  =  0.02  [cosh  120(zi  +  0.015)  cosh  120(x2  +  0.015)]"1  (180) 

The  domain’s  sides  were  0.1  meters  in  length.  The  unstructured  mesh  consisted  of 
700  nodes  and  1318  triangular  elements,  resulting  in  a  representative  element  size  of 
he  =  4.19  mm  (calculated  as  the  length  of  a  side  of  an  equilateral  triangle  with  an 
area  equal  to  the  mesh’s  mean  element  area).  The  same  physical  constants  used  in 
the  earlier  model  verification  were  retained. 

The  simulation  was  performed  using  the  Newmark  trapezoid  integration  method 
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(a)  Time  0 


(b)  Time  20 


(c)  Time  40 


(d)  Time  60 


(e)  Time  80 


Figure  40.  Time  slices  of  the  free  response  of  a  nonlinear  membrane. 
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and  the  Hermite  interpolation  method  at  time  steps  corresponding  to  CFL  num¬ 
bers  repeatedly  halved  from  0.8  to  0.025,  where  the  Courant-Friedrich-Levy  number 
CFL  =  cAt/he,  and  the  wave  speed  c  =  \J E / p  «  45  m/s.  Normalization  by  use 
of  the  CFL  number  is  for  intuitive  clarity  only,  since  both  methods  are  implicit  and 
therefore  technically  immune  to  C'FX-driven  stability  restrictions  [79].  The  duration 
of  the  simulation  was  80  time  steps  with  CFL  =  0.8.  For  convenience  when  compar¬ 
ing  time  series  with  different  time  increments,  we  will  refer  to  these  81  time  slices  as 
the  “reference  times.” 

The  remaining  discussion  will  focus  on  the  out-of-plane  response  u%  of  the  primary 
mesh  node  nearest  the  original  peak  at  (—0.015,  —0.015).  Since  an  analytical  solution 
or  a  proper  benchmark  solution  is  not  available,  a  baseline  was  derived  to  approximate 
the  true  solution.  Richardson  extrapolation  was  used  to  predict  the  solution  at  the 
reference  times  for  each  method  based  on  simulations  using  CFL  =  0.05  (1,280  time 
steps)  and  CFL  =  0.025  (2,560  time  steps).  The  nearly-coincident  extrapolated 
solutions  for  the  two  methods  are  depicted  in  Fig.  41.  They  are  essentially  converged, 
as  illustrated  by  plotting  the  difference  in  displacement  between  the  two  extrapolated 
solutions  in  Fig.  42.  The  largest  difference  was  on  the  order  of  1.0  x  10“6  m,  sufficiently 
small  when  considering  the  initial  peak’s  height  of  0.02  m.  To  avoid  biasing  the 
results  in  favor  of  either  technique,  the  baseline  for  the  following  discussion  was  then 
calculated  as  the  mean  of  the  two  method’s  extrapolated  solutions. 

As  expected  from  previous  results,  the  Hermite  method’s  solutions  remained  closer 
to  the  baseline  as  the  time  increment  was  increased.  Figure  43  compares  the  two 
methods  at  CFL  =  0.8.  The  dashed  line  of  the  Hermite  method  overlies  the  baseline 
for  much  of  the  simulation.  The  dispersion  of  the  Newmark  method  is  visible  as  the 
solution  progressively  lags  the  baseline.  A  convergence  study  was  performed  using 
the  L2  norms  of  the  errors  at  the  reference  times,  and  the  results  are  plotted  in 
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Figure  41.  Converged  extrapolated  solutions  for  the  Hermite  and  Newmark  trapezoidal 
methods. 


Figure  42.  Difference  in  displacement  between  the  extrapolated  solutions  of  the  Her¬ 
mite  and  Newmark  trapezoidal  methods. 
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Fig.  44.  The  observed  convergence  rates  of  the  two  methods  between  the  two  finest 
time  meshes  were  2.06  for  the  Hcrmite  method  and  1.99  for  the  Newmark  method. 
Thus  the  results  of  this  case  study  corroborate  those  of  the  verification  procedure: 
the  observed  rates  of  convergence  of  both  methods  are  two,  but  the  Hermite  method 
produces  significantly  smaller  errors. 


Figure  43.  Comparison  of  the  model  solutions  at  CFL  =  0.8  to  the  baseline  solution. 


Figure  44.  Convergence  plots  of  the  dynamic  nonlinear  membrane. 

At  the  time  step  sizes  presented  here,  the  conditional  stability  of  the  Hcrmite 
method  was  not  a  factor.  Performance  of  both  time-marching  methods  began  to 


123 


suffer  significantly  at  larger  time  increments  -  the  cost  of  each  time  step  increased 
dramatically  while  accuracy  deteriorated.  The  increased  cost  was  attributed  to  the 
decreased  accuracy  of  the  initial  guesses.  The  nonlinear  solver  required  more  itera¬ 
tions  to  converge,  and  convergence  became  less  assured.  In  onr  experience  the  best 
combination  of  computational  speed  and  solution  accuracy  was  achieved  using  the 
Hcrmite  method  with  smaller  time  steps,  where  the  fixed-point  iterations  converged 
rapidly.  A  more  rigorous  cost  comparison  would  be  required  to  quantify  the  differ¬ 
ences  between  the  two  methods  and  determine  optimal  time  step  sizes. 

Summary. 

This  section  presented  the  performance  of  a  dynamic  model  for  geometrically 
nonlinear  membranes.  The  model  incorporated  the  two  previously  proposed  dis¬ 
cretization  methods:  a  jerk  constraint-based  Hermite  time  interpolation  method,  and 
a  staggered-grid  point  collocation  method  with  grouped  nonlinear  terms.  The  Her¬ 
mite  method  used  the  physical  meaning  of  the  rate  of  change  of  acceleration  of  the 
structure  to  define  the  jerk  constraints.  The  Newmark  trapezoid  method  was  selected 
for  comparison  because  it  is  well  understood  and  is  the  generally  accepted  standard 
for  structural  engineering.  The  hybrid  method  was  evaluated  during  verification  and 
found  to  fill  the  middle  ground  between  the  Newmark  and  Hermite  methods  -  it 
enabled  larger  time  step  sizes  than  the  pure  Hermite  method,  and  the  accuracy  for 
out-of-plane  displacements  was  better  than  the  pure  Newmark  scheme.  Two  evalua¬ 
tion  cases  were  performed,  one  forced  and  one  free,  by  changing  time  increment  sizes 
for  a  given  mesh.  First,  the  Method  of  Manufactured  Solutions  was  employed  on  a 
hexagonal  domain  to  verify  the  code  and  calculate  rates  of  convergence.  Second,  a 
free-response  case  study  was  evaluated  on  a  square  domain.  The  nearly-converged  ex¬ 
trapolated  solutions  of  the  two  time-marching  methods  were  used  to  create  a  baseline 
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solution.  The  results  generally  confirmed  the  behavior  of  the  Hermite  time  method: 
minimal  frequency  distortion  and  improved  accuracy  versus  the  Newmark  method. 
In  a  notable  exception,  the  observed  rate  of  convergence  for  the  Hermite  method  fell 
from  four  in  previous  evaluations  to  two  for  this  membrane  model.  This  rate  was 
observed  across  both  evaluation  cases  and  despite  different  norms.  Overall,  for  this 
particular  application,  the  combination  of  the  two  proposed  discretization  schemes 
was  feasible  and  suitable  for  obtaining  significantly-improved  accuracy. 
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V.  Conclusions 


The  objective  of  this  research  was  to  devise  and  evaluate  novel  discretization 
schemes  for  the  dynamic  simulation  of  membranes.  The  study  was  divided  into  two 
primary  focus  areas:  (1)  linear  membranes,  and  (2)  geometrically  nonlinear  mem¬ 
branes,  as  defined  in  Chapter  II.  Different  techniques  were  developed  and  tested  for 
the  two  cases.  Throughout  this  study,  the  Newmark  trapezoid  method  was  selected 
for  comparison  because  it  is  well  understood  and  is  the  generally  accepted  standard 
for  structural  engineering.  In  this  final  chapter,  the  research  will  be  summarized 
and  discussed,  contributions  will  be  stated,  and  avenues  for  future  research  will  be 
suggested. 

Linear  Membrane  Dynamics 

The  Simultaneous  Time-Continuous  Galerkin  (STCG)  method  was  applied  to  a 
classical  linear  membrane  model.  The  method  featured  trilinear  space-time  elements 
and  a  mixed  formulation.  The  entire  space-time  domain  was  discretized,  and  the  sim¬ 
ulation’s  entire  solution  was  obtained  by  solving  the  single  large,  sparse,  linear  system. 
Numerical  studies  indicated  observed  second-order  rates  of  convergence  in  both  space 
and  time,  which  is  typical  of  linear  elements.  Bounded  numerical  instabilities  in  the 
form  of  high-frequency  oscillations  were  present.  However,  the  lower  modes  were  ac¬ 
curately  represented,  so  a  post-solution  smoothing  procedure  was  demonstrated  that 
effectively  eliminated  high-frequency  oscillations.  Unique  to  this  method,  when  the 
time  steps  exceeded  a  critical  frequency,  the  solution  rapidly  collapsed  to  its  resting 
configuration. 

The  spurious  modes  were  likely  caused  by  the  combination  of  interpolations  used 
for  discretizing  the  mixed  form.  As  discussed  in  [146],  mixed  forms  require  appropriate 
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matching  of  interpolations  to  satisfy  stability  conditions.  Fortunately,  for  the  scheme 
presented  here,  the  spurious  oscillations  were  bounded  and  not  destructive  to  the 
lower-frequency  modes.  Post-solution  filtering  removed  much  of  the  unwanted  content 
with  slight  dissipation  of  the  peaks  of  the  lower  modes. 

Nonlinear  Membrane  Dynamics 

The  dynamic  nonlinear  membrane  model  was  composed  of  two  independently  de¬ 
veloped  schemes:  a  jerk  constraint-based  Hermite  time  interpolation  method  detailed 
in  Section  IV,  and  a  staggered-grid  point  collocation  method  with  grouped  nonlin¬ 
ear  terms  described  in  Section  IV.  Each  scheme  was  rigorously  verified  and  validated 
through  the  use  of  cases  with  known  analytical  solutions,  manufactured  solutions,  and 
experimental  data.  Finally,  the  two  schemes  were  combined  to  form  a  pure  collocation 
model,  which  was  analyzed  in  Section  IV. 

The  Hermite  method  used  the  physical  meaning  of  the  rate  of  change  of  accel¬ 
eration  to  define  the  jerk  constraints  and  the  resulting  unique  quintic  polynomial 
interpolation  of  displacement  during  a  time  step.  The  single-step,  implicit  scheme 
was  shown  to  be  significantly  more  accurate  than  the  most  common  second-order 
Newmark  schemes  (central  difference  and  trapezoid).  It  was  proved  that  the  New- 
mark  scheme  is  a  subset  of  the  proposed  scheme,  and  is  obtained  by  appropriately 
defining  the  jerk  constraints.  It  was  also  proved  that  the  scheme  is  symplectic  for  a 
simple  harmonic  oscillator.  Period  elongation  was  derived  analytically  for  a  simple 
harmonic  oscillator  and  shown  to  be  minimal.  Through  a  classical  linear  analysis,  the 
scheme  was  found  to  be  conditionally  stable,  and  the  limits  of  stability  were  analyti¬ 
cally  derived.  The  regions  of  instability  were  found  to  be  above  the  Nyquist  frequency. 
The  method  is  not  suitable  for  mathematically  stiff  systems.  However,  a  novel  error 
estimation  scheme  was  developed  which,  when  used  to  control  time  step  sizes,  enabled 
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the  method  to  successfully  solve  a  stiff  system  in  a  numerical  demonstration. 

The  Hermite  time  interpolation  method  shows  promise  for  general  application  to 
differential  equations,  as  shown  by  the  variety  of  numerical  examples.  The  method 
proposed  in  this  paper  provides  a  unique  avenue  of  investigation.  Since  jerk  is  not 
continuous  across  time  steps,  the  jerk  constraints  can  be  modified  or  optimized  sepa¬ 
rately  for  each  time  step  or  degree  of  freedom  to  achieve  other  desired  effects  such  as 
algorithmic  damping  or  explicit  energy  conservation.  This  feature  was  illustrated  by 
the  hybrid  method  of  Section  IV,  where  the  jerk  constraints  were  defined  differently 
in-plane  versus  out-of-plane. 

The  proposed  spatial  discretization  scheme  accounts  for  both  in-plane  and  out- 
of-plane  displacements  without  carrying  rotational  degrees  of  freedom.  Based  on 
the  Method  of  Weighted  Residuals  point  collocation  approach,  it  employs  a  staggered 
grid,  grouping  of  nonlinear  terms,  and  polygon  shape  functions  in  a  strong-form  point 
collocation  formulation.  Additional  demonstrated  capabilities  include  varying  mem¬ 
brane  thickness  and  follower  forces,  as  demonstrated  in  Section  IV.  The  observed  rate 
of  convergence  was  two  for  both  displacement  and  strain.  Validation  against  exper¬ 
imental  data  in  the  literature  showed  the  method  to  be  accurate  until  hyperelastic 
effects  begin  to  appear. 

The  point  collocation  scheme  is  simpler  to  formulate  than  conventional  nonlin¬ 
ear  finite  element  approaches.  Element  integration  is  avoided  entirely,  a  significant 
benefit  since  with  classical  Newton  iterations,  each  integration  point  (with  at  least 
three  Gaussian  integration  points  per  triangular  element)  must  be  computed  for  each 
iteration  of  each  time  step.  The  group  formulation  permits  the  same  treatment  in 
all  three  axes,  again  simplifying  coding.  Overall,  the  framework  of  the  approach  is 
highly  modular  and  flexible.  Any  given  step  can  be  performed  by  interchangeable 
subroutines.  The  residual  subroutine  readily  accepts  different  strain-displacement  or 
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material  models  (including  nonlinear  models). 

Lastly,  the  Hermite  time  interpolation  method  and  the  point  collocation  spatial 
scheme  were  combined  to  create  a  novel  dynamic  membrane  model  based  purely  on 
collocation  techniques.  Two  evaluation  cases  were  considered,  one  forced  and  one 
free.  Convergence  studies  were  performed  to  observe  the  convergence  rates  with  re¬ 
spect  to  time  step  size.  First,  the  Method  of  Manufactured  Solutions  was  employed 
on  a  hexagonal  domain  to  verify  the  code  and  calculate  rates  of  convergence.  Sec¬ 
ond,  a  free-response  case  study  was  performed  to  evaluate  the  model’s  ability  to 
capture  traveling  wave  phenomena.  The  converged  extrapolated  solutions  of  the  two 
time-marching  methods  (Hermite  and  Newmark  trapezoid  methods)  were  averaged  to 
create  a  baseline  solution.  Time  step  sizes  were  then  increased  to  observe  how  rapidly 
the  solutions  from  the  two  methods  deviated  from  the  baseline  solution.  The  results 
generally  confirmed  the  behavior  of  the  Hermite  time  method  observed  in  Section  IV: 
minimal  frequency  distortion  and  improved  accuracy  versus  the  Newmark  method. 
It  was  notable  that  the  observed  rate  of  convergence  for  the  Hermite  method  fell 
from  four  in  the  evaluations  discussed  in  Section  IV  to  two  for  the  dynamic  mem¬ 
brane  simulation  discussed  in  Section  IV.  Overall,  for  this  particular  application,  the 
combination  of  the  two  proposed  discretization  schemes  was  feasible  and  suitable  for 
obtaining  significantly-improved  accuracy. 

Contributions 

This  research  effort  involved  the  assessment  of  a  wide  variety  of  numerical  meth¬ 
ods  for  use  in  membrane  simulation,  and  ultimately  culminated  in  a  novel,  fully 
three-dimensional,  highly  modular  membrane  model  that  accurately  and  efficiently 
simulates  the  dynamics  of  geometrically  nonlinear  membranes.  It  extends,  rigorously 
investigates,  and  incorporates  specific  techniques  that  are  not  typically  used  in  struc- 
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tural  dynamics.  The  process  of  developing  and  evaluating  this  model  has  led  to  the 
following  specific  contributions: 

1.  Hcrmite  time  interpolation  method 

•  First  formulation  by  strongly  enforcing  jerk  constraints.  This  approach 
offers  the  potential  for  future  optimization  by  the  selection  of  different 
constraint  formulas. 

•  First  direct  link  proven  between  the  jerk  constraint-based  Hermite  inter¬ 
polation  method  and  the  Newmark  family  of  methods. 

•  First  precise  analytical  expression  for  the  linear  stability  limits. 

•  First  precise  analytical  expression  for  dissipation  of  a  simple  harmonic 
oscillator. 

•  First  local  error  estimation  technique  for  a  Hcrmite  time  interpolation.  The 
accurate  local  error  estimation  capitalized  on  the  intra-element  solution 
accuracy  by  using  mid-time-step  solutions  to  estimate  error. 

•  First  solution  of  a  stiff  system.  The  local  error  estimation  technique  ef¬ 
fectively  controlled  time  step  size  to  accurately  solve  a  numerically  stiff 
system. 

2.  Point  collocation  model 

•  First  derivation  and  application  of  spatial  derivatives  of  the  polygon  shape 
functions  proposed  in  [29]. 

•  First  use  of  the  group  finite  element  method  in  a  point  collocation  model. 

•  First  use  of  a  staggered  two-dimensional  grid  approach  for  nonlinear  PDE 
solution. 
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3.  Geometrically  nonlinear  membrane  model 


•  First  to  incorporate  Hcrmite  time  interpolation 

•  First  to  incorporate  a  staggered-grid  point  collocation  spatial  discretization 

Each  of  the  proposed  discretization  schemes  stands  on  its  own  as  an  effective 
numerical  method  for  solving  nonlinear  ordinary  or  partial  differential  equations. 
Therefore,  the  contributions  of  this  study  extend  well  beyond  Micro  Air  Vehicles, 
membranes,  and  structural  analysis. 

Recommendations  for  Future  Work 

Further  investigation  of  the  STCG  concept  should  start  with  determining  the  pre¬ 
cise  source  of  the  high-frequency  noise  and  mitigating  it  intrinsically,  so  post-solution 
filtering  is  not  necessary.  As  discussed  in  Section  III,  the  results  of  the  static  verifi¬ 
cation  indicated  a  mild  instability  in  the  mixed-form  solution.  Using  a  displacement 
formulation  rather  than  a  mixed  formulation  in  space  could  test  this  source  of  oscilla¬ 
tions,  as  well  as  reduce  the  number  of  degrees  of  freedom  for  the  global  system.  After 
addressing  the  high-frequency  noise,  the  temporal  stability  behavior  should  be  ana¬ 
lytically  derived.  In  particular,  eigenvalue  analysis  may  explain  the  tendency  of  the 
solution  to  collapse  when  the  time  step  size  exceeds  a  certain  threshold.  Finally,  the 
computational  efficiency  could  be  compared  to  reference  schemes  such  as  the  linear 
Newmark  and  explicit  Runge  Kutta  methods. 

Future  studies  of  the  Hermite  interpolation  method  can  address  the  following  three 
issues:  (1)  stability,  (2)  numerical  damping,  and  (3)  computational  efficiency.  Since 
only  linear  stability  has  been  addressed  in  this  research,  nonlinear  stability  behavior 
could  be  characterized.  The  conditional  stability  of  the  method  is  a  significant  draw¬ 
back  for  some  applications.  Efforts  to  expand  the  stability  regimes,  both  intrinsically 
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and  using  feedback  controls,  could  increase  the  utility  of  the  method.  A  component 
of  stability  is  the  numerical  damping  of  high-frequency  content.  The  excellent  energy 
preservation  is  a  useful  feature  of  the  method;  however,  the  use  of  a  non-dissipative 
algorithm  for  a  hyperbolic  problem  raises  the  spectre  of  error  propagation  and  build¬ 
up,  which  can  cause  a  solution  to  fail.  There  will  likely  arise  a  trade-off  between 
accuracy  and  damping  effectiveness,  though  by  the  time  the  damping  is  needed,  the 
local  solution  may  have  already  departed  significantly  from  the  exact  solution. 

For  the  point  collocation  spatial  discretization  scheme  of  Section  IV,  further  ef¬ 
forts  could  investigate  the  sensitivity  of  the  solution  to  mesh  geometry,  in  particular 
through  analytical  a  priori  error  estimates  and  determination  of  the  scheme’s  order  of 
accuracy.  Also,  the  order  of  the  method  could  be  derived  analytically  so  the  a  priori 
estimate  could  be  verified  by  the  numerical  convergence  studies.  Free  edge  boundary 
conditions  for  the  strong  form  formulation  could  potentially  be  difficult  to  formulate, 
so  a  parallel  development  of  the  overall  concept  (staggered  grid  with  a  group  FE 
formulation)  from  the  weak  form  would  probably  illuminate  a  practical  path  of  pur¬ 
suit.  The  computational  expense  of  using  the  proposed  scheme  with  a  Jacobian-free 
nonlinear  solver  could  be  compared  to  a  classical  linearized  finite  element  formulation 
using  Newton  iterations.  Finally,  application  of  the  method  to  other  systems  such 
as  hyperelastic  membranes  and  thin  plates  offers  many  interesting  opportunities  for 
future  work. 
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